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I .  INTRODUCTION 


The  objective  of  this  analysis  has  been  to  develop  a  finite  ele¬ 
ment  model  to  analyze  dynamically  loaded,  laminated,  orthotropic  thick 
plates  with  efficient  numerical  methods  of  solution.  It  is  intended 
to  apply  the  analysis  to  relatively  large  magnitude  forces  which  are 
rapidly  changing  with  time,  consequently  the  materials  are  assumed  to 
be  represented  by  elastic-plastic  and  elastic-viscoplastic  models. 

The  finite  element  model  assumes  large  transverse  deflections  of  the 
plate  which  lead  to  nonlinear  strain- displacement  relations.  Because 
the  analysis  is  intended  for  predicting  inter-laminar  stresses  and 
thickness  changes,  basic  plate  theory  will  not  suffice.  However,  full 
three-dimensional  analysis  is  much  too  cumbersome  to  use,  so  this  me¬ 
thod  makes  use  of  the  plate-like  characteristics,  but  still  allows  for 
three-dimensional  response.  The  finite  element  model  leads  to  a  set 
of  dynamic  matrix  equations  representing  the  nodal  plate  displacements. 
The  basic  approach  to  the  solution  of  these  equations  consists  of  an 
incremental  approach  in  the  time  domain.1 

The  method  of  attack  is  to  first  calculate  the  internal  forces 
from  the  incremental  stresses  and  deflections  of  the  previous  time  in¬ 
crement  using  the  stiffness  matrix.  The  external  load  is  input  and, 
using  the  theory  of  virtual  work,  it  is  transformed  to  external  nodal 
forces.  Using  Newmark's  beta  finite  difference  technique,  the  de¬ 
flections  for  the  new  time  increment  are  calculated.  From  these  new 
deflections  the  stresses  are  calculated  and  these  stresses  are  checked 
to  see  if  elastic-plastic  or  elastic-viscoplastic  yielding  has  occurred. 
The  time  is  then  incremented  up  and  this  technique  marches  on  for  the 
desired  length  of  time.  Figure  1  shows  a  flow  diagram  of  the  computer 
program  with  appropriate  section  references. 

This  report  will  first  review  the  basic  finite  element  formulation 
of  the  problem  which  will  include  the  development  of  the  stiffness  ma¬ 
trix  and  the  stress  calculations.  It  will  then  discuss  the  development 
of  the  external  nodal  forces,  the  dynamic  analysis  with  a  variable  time 
increment,  and  the  orthotropic,  elastic-plastic  and  elastic-viscoplastic 
analysis.  The  remainder  of  the  analysis  will  be  devoted  to  the  methods 
used  to  increase  the  size  of  the  time  increment  in  the  time  integration 
and  the  element  equilibrium  equations  that  resulted  from  this  increase 
in  time  increment  size.  Finally,  a  numerical  example  will  be  discussed, 
and  the  accuracy  of  this  method  of  analysis  will  be  compared  to  actual 
experimental  results. 
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Review  of  Finite  Element  Formulation  for  Laminated  Plates 


Because  of  the  physical  characteristics  of  a  thick  laminated  plate, 
it  is  natural  to  use  Cartesian  coordinates  and  to  align  the  plane  of 
the  plate  with  the  xi  -  X2  coordinate  plane  and  the  thickness  with  the 
X3  coordinate  axis.  The  shape  of  the  finite  element  will  then  be  a  gen¬ 
eral  quadrilateral  defined  in  the  xi  -  X2  coordinate  plane  with  a  con¬ 
stant  thickness  in  the  X3  coordinate  axis.  To  insure  a  homogeneous 
element,  there  should  be  as  many  elements  stacked  in  the  X3  direction  as 
there  are  different  laminated  layers.  A  typical  numbering  system  and 
axis  location  is  shown  in  Figure  2.  The  mass  of  each  element  is  lumped 
at  the  nodes.  The  actual  calculation  of  the  mass  at  each  node  is  cal¬ 
culated  by  bisecting  each  side  of  the  triangular  element  and  joining 
these  center  points  to  the  center  node  which  defines  the  apex  of  the  four 
triangular  elements  in  Figure  3.  Thus,  the  original  quadrilateral  has 
been  divided  into  four  smaller  ones  which  contain  two  of  the  eight  nodes 
of  the  large  element.  The  mass  of  each  smaller  quadrilateral  portion  of 
the  element  is  found  by  multiplying  the  thickness  of  the  element  by  the 
density  by  the  area  of  the  smaller  quadrilateral  element,  and  the  result 
is  distributed  equally  to  the  two  nodes.  The  mass  matrix  is  obtained  by 
summing  masses  at  each  structural  node  from  the  adjacent  elements.  The 
resulting  matrix  is  diagonal. 

For  ease  in  calculating  the  stiffness  matrix,  each  quadrilateral 
element  is  subdivided  into  four  triangular  elements  with  the  center  node 
not  carrying  any  load  so  that  it  can  later  be  removed  by  static  conden¬ 
sation  when  the  triangular  elements  are  summed  to  form  the  quadrilateral 
element.  A  typical  element  and  sub-element  nodal  numbering  system  is 
shown  in  Figure  3. 


Displacement  Functions 

The  first  step  in  finite  element  formulation  is  to  choose  a  dis¬ 
placement  function.  Since  each  traingular  element  has  three  degrees  of 
freedom  at  each  node,  the  assumed  displacement  function  for  each  ele¬ 
ment  should  have  eighteen  generalized  coordinates.  The  displacement 
function,  chosen  to  be  linear  in  the  x^  -  X2  plane  and  to  vary  linearly 
in  the  X3  direction,  is  shown  below: 
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u3  -  a7  Vl  +  V2  +  X3  (B7+B8X1+B9X2) 

where  u^  (i=l,3)  are  the  Cartesian  displacement  components  in  terms  of 
x^  directions  and  otj  and  (3j  (j=l,9),  the  eighteen  unknown  generalized 
coordinates.  Writing  Equations  1.1  at  each  of  the  six  nodes  of  a  trian¬ 
gular  element  results  in  eighteen  equations  which  can  be  written  as  the 
matrix  shown  below: 
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where  the  forward  superscripts  on  ui  and  quantities  refer  to  one  of 
the  six  nodes  and  t  represents  the  element  thickness.  Equation  1.2  can 
be  written  in  a  short  form  as: 

{6}  =  [E]  {ot}  1.3 

Solving  for  the  global  coordinates  gives: 

{a}  =  [E]_1{6>  1.4 

Because  there  are  so  many  zeros,  matric  [E]  can  be  inverted  by  hand  in 
the  following  manner.  First,  divide  each  column  matrix  in  half  so  there 
is  an  upper  and  lower  generalized  coordinate  matrix  (a, 3)  and  an  upper 
and  lower  nodal  displacement  matrix  (6  ,  6fl).  This  changes  Equation  1.3 
to:  u  1 
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[e2] t[E2] 

{3/ 

V  J 

where  the  [E]  matrix  has  been  divided  into  four  parts.  Inverting  Equa¬ 
tion  1.5  by  parts  yields: 

{a}  =  [E1 ] -1{6  }  1.6 

L  1J  u 

{3>  =  i[E2J-1{{6Jl}  -  [E2]{a}}  1.7 

Inserting  Equation  1.6  into  1.7  yields: 
tB}  - 

Using  Equations  1.6  and  1.8,  it  is  apparent  that  Equation  1.4  can  be 
written  as : 
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Using  the  following  notation: 
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2k^  ,  =  det 

t  or  b 


i  1  1 

1  Xj  x2 

1  3Xj  jx2 

_  m  m 
1  X1  x2 


where  i,j,m  are  cyclic  permutations  either  of  1  to  3  for  b  (bottom  of 
the  element)  or  4  to  6  for  t  (top  of  the  element),  the  individual  matri¬ 
ces  are  as  follows: 
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12  3 
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0  a^  0  0  a2  0  0  a^  0 
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a.  0  0  a,.  0  0  a.  0  0 

4  5  6 
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Strain-Displacement  Relation 

The  next  step  in  the  finite  element  formulation  is  to  calculate  the 
internal  work  due  to  a  virtual  change  in  nodal  displacement.  The  strains 
must  first  be  found  in  terms  of  the  generalized  coordinates.  The  non¬ 
linear  strain  displacement  relations  in  tensor  notation  are: 


'ij 


2l  ij  j,: 


1.13 


Making  use  of  the  fact  that  this  analysis  is  for  a  plate,  it  is  assumed 
that  deflections  and  rotations  out  of  the  plane  of  the  plate  are  large 
compared  to  those  in  the  plane  of  the  plate,  or  symbolically: 
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1.14 

Ul,3 

,U2,3  >  U3,3 

Using  the  order  of  magnitude  defined  in  Equations  1.14  reduces  Equations 
1.13  to  the  following  form: 
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Since  there  are  linear  and  nonlinear  terms.  Equation  1.15  can  be  re¬ 
written  in  matrix  notation  as : 
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^en£^  2^U3,1,U3,2,U1,3  +  U2,3’U3,1U3,2,0’°^  1,18 

Performing  the  indicated  operations  in  Equation  1.17  on  1.1  results  in 
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In  matrix  notation  Equation  1.19  reduces  to: 

{e^}  =  [Q]  {a}  1.20 

Substituting  Equation  1.4  into  1.20  results  in: 

{e£>  =  [Q]  [E] _1{6}  1.21 

Taking  virtual  changes ,  this  can  be  written  as : 

d  {e£>  =  [B£]  d  {6}  1.22 

where: 

[B£]  =  [Q]  [E ] _ 1  1.23 

and  the  symbol  d  in  front  of  a  matrix  represents  a  virtual  change. 
Disregarding  the  zero  terms.  Equation  1.18  can  be  written  as: 


11 


r  % 

Elln£, 

U3,I  ° 

0 

0 

U3,l 

.  Z22nl 

1 

2 

0  u3,2 

0 

0 

- 

U3,2  ■ 

e33nA 

0  0 

Ul,3 

U2,3 

Ul,3 

i  £12nAj 

_U3,2  U3.I 

0 

0 

.  U2 , 3  J 

In  matrix  notation  Equation  1.24  reduces  to: 

i[A]  {0}  1.25 

Since  the  virtual  changes  are  considered  in  the  finite  element  method, 
it  is  necessary  to  investigate  virtual  changes  in  the  nonlinear  strain 
which  can  be  written  as  follows: 

d  {en!}  =  \  [A]  d  {0}  +  |d  [A]  {0}  1.26 

Clearly  from  Equation  1.24: 

[A]  d  {0}  =  d  [A]  {0}  1.27 

From  Equations  1.26  and  1.27,  the  following  result  is  true: 

d  {en!}  =  [A]  d  {0}  1.28 

Performing  the  indicated  operations  in  Equation  1.24  on  1.1  results  in 
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Equation  1.31  can  be  written  in  a  short  form  as  follows: 

d  {en£*  =  [Z]  d  ®  1,32 

where  the  matrices  [a],  [Z]  and  d{a}  are  symbolic  representation__of  the 
three  matrices  on  the  right  hand  side  of  Equation  1.31.  Since  [a]  are 
constants  which  can  be  written  in  terms  of  the  rows  of  [E]-*  correspon¬ 
ding  to  the  correct  aj  and  {6}  given  in  Equation  1.4,  Equation  1.32  can 
be  rewritten  as: 

d  *enp  =  ®  [Z]  [E]  d  {6}  1.33 

where  [E]  also  consists  of  the  rows  corresponding  to  the  correct  a.  in 
{a}.  This  is  written  as  the  normal  strain-displacement  relationship 
below: 

d  {En*}  ■  tB„d  d  t6>  !-34 

where: 

[Bn£]  =  [a]  [Z]  [E]  1.35 

Note  that  [B  ]  is  dependent  on  displacements  from  [  ] .  From  Equations 
1.16,  1.22  and  1.34,  the  complete  strain-displacement  relationship  be¬ 
comes  : 


d  {e}  =  [B£]  d  {&}  +  [B  £]  d  {6}  1.36 

which  in  general  is  written  as : 

d  {e}  =  [B]  d  {6}  1.37 

where : 

[B]  =  [B £]  +  [BnA]  1.38 


Stress-Strain  Relationship 


The  only  other  calculation  to  be  made  before  the  virtual  work  is 
evaluated  is  the  stress-strain  relationship.  Since  this  analysis  deals 
with  orthotropic  material,  only  the  orthotropic  relationship  will  be 
discussed  since  in  the  limit  these  equations  become  isotropic  when  the 
modulus  of  elasticity  and  Poisson* s  ratio  are  input  into  the  following 
equations : 
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en*  =  g“  [Pu  '-v^Paa’-ViaPaa'  ] 

£22*  =  ==■  [P2  2  *  ”^2  lPl 1  * “^2  3P3  3  *  3 

C2 


E3  3  ' 


“V3  1 CJ 1  1  ’  -V32CJ22  *  ] 


1.39 


El2 


f 


Pi  2' 

2G12 


^2  3  ' 


P2  3* 

2G2  3 


Els' 


Pis' 

2Gi  3 


where  is  the  modulus  of  elasticity  of  the  x^  direction,  v—  and 
G^j  are  the  Poisson's  ratio  and  shear  modulus  along  the  x^-xj  plane. 

ItJ should  be  noted  that  the  prime  in  these  equations  refers  to  the  local 
numbering  system.  Due  to  the  plane  of  symmetry  for  orthotropic  mater¬ 
ials,  the  following  is  true: 


v12  _ 
Ei 


Vj  3  _  V3 1 

Ei  E3 

1.40 

V2 3  _  V32 

E2  "  E3 

Letting: 

L,  =  I-V32V2  3-V12V2 1 -Vi  2V2  3V3 1 -Vi  3V32V21  1.41 


the  inverse  of  Equation  1.39  is: 
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al 1  f “  (1-V3 2V2 3 ) E 1  £ 1 1  f  +  (Vi 2+Vi 3V2  3) E 2^2 2 ,  +  (Vi 3  +  V12^2  3)E3e3  3  f ] 


a2  2  1  -  “■[  (V2 i+V2  3^3l)El£li!  +  (l“Vi3V3i)E2£2  2f  +  (V2  3+^2  lVi3)E3£33f] 


a3  3  f“  £*[  (V3  1+V3  2V2  l)Ei£i  1  f  +  (V  3  2+V  3  1 V 1  2)  E2  £2  2  1  +  (1  "v  1  2  v2  1 )  E  3  £  3  3  f  ] 


Cf  1  2  * 

2G12£i2f 

1.42 

a2  3  * 

=  2G23£23f 

cri3f 

=  2Gi 3£i 3  f 

In  matrix  notation  Equation 

1.42  is: 

{o'} 

=  [C']  {£'} 

1.43 

Often,  though,  the  axes  for  which  the  material  properties  are  defined  do 
not  align  with  those  of  the  global  coordinate  system.  When  this  happens 
the  stresses  of  interest  are  those  of  the  global  system.  Letting  the 
global  system  have  the  unprimed  stresses  and  strains,  an  equation  simi¬ 
lar  to  Equation  1.43  can  be  written  as  follows: 

{a} 

=  [C]  {£} 

1.44 

In  order  to  transform  {o'}  to  {a}  the  direction  cosines  between  x^  and 

x^  can  be  used.  We  define  directional  cosines  by  a • •  as  follows: 

1 


a.  .  =  cos  (x. ,x1)  1 . 45 

l]  1  J 


and  therefore  from  stress  transformation  relations,  it  follows  that: 


V  ■ 

and: 

£,  o '  =  a.,  a.  „e.  . 

k£  lk  ij 


1.46 


1.47 
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Equations  1.46  and  1.47  can  be  rewritten  in  matrix  form  as: 

{a'}  =  [R]  {a}  1.48 

and : 

{£'}  =  [R]  {el  1.49 

where : 


an 

ah 

all 

2ana2i 

2a3  xa2 1 

2a  1 1«3 1 

a?  2 

a!  2 

a  3  2 

2aj2a22 

2a32a22 

2a12a32 

a?  3 

al3 

a  3  3 

2a13a23 

2a3  3a2  3 

2oti  3a3  3 

al lal 2 

a2 la22 

a3 !a32 

ai  ia2  2+a2  iai  2 

ai ias2+a3  jaj  2 

a2 ia32+a3  ja22 

ai 2«i 3 

&2  2^2  3 

a32a3  3 

al  2a2  3  +  a2  2ai  3 

a22«3  3  +  OC3  2a2  3 

ai 2a3  3+a32ai 3 

ana13 

a2ia23 

a3ia33 

aua23+a2iai3 

a2  ia3  3+a3  ^2  3 

anaja+agian 

The  work  done  in  either  coordinate  system  must  be  equal,  therefore: 

d  {e}T{a}  =  d  {e'}T  {a'}  1.51 

Substituting  Equations  1.43,  1.44  and  1.49  into  1.51  results  in: 

d  {e}T  [C]  {e}  =  {e}T  [R]T  [C']  [R]  {e}  1.52 

Therefore,  from  Equation  1.52  the  stress-strain  matrix  in  global  coor¬ 
dinates  can  be  written  in  terms  of  the  material  stress-strain  relation¬ 
ship  and  the  transformation  matrices.  This  can  be  written  as: 

[C]  =  [R]T  [C']  [R]  1.53 

In  the  computer  program  the  transformation  matrix  is  found  by  two 
separate  transformations.  Figure  4  shows  the  orientation  of  the  coordi¬ 
nate  systems  and  the  two  angles  needed  to  describe  the  transformation. 


18 


Note  that  lies  in  the  Xj-X2  plane.  Using  the  a  and  3  described 
here,  not  to  be  confused  with  the  subscripted  generalized  coordinates, 
the  transformation  matrices  can  be  calculated  by  the  appropriate  in¬ 
sertion  into  the  direction  cosines. 


Internal  Force 


It  is  now  possible  to  calculate  the  virtual  change  in  the  internal 
work  of  the  structure  due  to  a  virtual  change  in  the  nodal  displace¬ 
ments.  This  can  be  written  as: 


d  W 


d  {e}T  {a}  dV 
v 


1.54 


Using  Equations  1.37 


in  1.54,  the  internal  work  is  written  as: 


d  W 


d  {6} 


[B]  {a}  dV 


1.55 


In  addition  to  the  work  done  by  internal  work  of  the  finite  element, 
there  is  external  work  done  by  the  forces  at  the  eight  corners.  De¬ 
noting  these  forces  by  the  matrix  {f},  the  external  virtual  work  is: 


d  W  =  d  {6}T  {f} 
E 


1.56 


Since  external  and  internal  forces  on  each  element  are  in  equilibrium, 
it  follows  from  virtual  work  that: 


d  WT  +  d  W„  =  0 

I  E 

Therefore,  from  Equations  1.55,  1.56  and  1.57: 


1.57 


{f} 


f  T 

[B]  {a}  dV 

'  V 


1.58 


The  forces  acting  on  any  node  resulting  from  the  adjacent  elements 
are  obtained  from  Equation  1.58.  This  summation  will  result  in  a  ma¬ 
trix  of  internal  forces  denoted  by  {Fj}. 


The  next  step  in  the  analysis  is  to  obtain  the  nodal  equilibrium 
equations  for  the  total  structure.  Because  this  is  a  dynamic  problem 
the  inertial  effects  must  be  considered.  The  equation  of  motion  is: 
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1.59 


[M]  {A}  =  {Fj}  +  {Fe> 


where  [M]  is  the  diagonal  mass  matrix  previously  described,  {A}  is  the 
global  displacement  acceleration,  and  {Fg}  is  the  concentrated  nodal 
force  matrix  to  be  described  in  Section  2.  The  solution  of  Equation 
1.59  will  be  discussed  in  Section  3. 
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II.  EXTERNAL  FORCE 


Because  of  the  varying  types  of  distributed  loads  that  could  be 
used  in  a  problem  of  this  type,  a  technique  was  developed  that  only 
needed  the  input  of  the  distributed  external  load  at  the  nodal  locations 
that  the  load  acts  on.  This  is  normally  accomplished  by  adding  a  sub¬ 
routine  to  the  program  which  gives  that  data.  Using  the  theory  of  vir¬ 
tual  work  it  is  possible  to  transform  these  external  distributed  loads 
to  the  effective  concentrated  forces  {Fg}  which  are  used  in  Equation 


1.59. 


Since  the  top  of  the  plate  which  carries  the  distributed  loads  is 
composed  of  triangular  elements,  as  shown  in  Figure  3,  the  simplest  re¬ 
presentation  of  the  distributed  load  L  over  the  plate  is : 


L  =  £1  +  £2x1  +  £3x2 


2.1 


where  £.  are  some  suitable  parameters  defining  the  linear  distribution. 
Writing  Equation  2.1  at  the  three  nodes  results  in: 


\ 

L1 

r-H 

1  y 

X1 

1x2 

t  " 

■L2 

» 

1 

2y 

X1 

2x2 

* 

£2 

vL3. 

1 

3  v 

X1 

3X2 

£3 

where  forward  superscripts  in  x.  refer 
notation  defined  by  Equation  l.lo,  the 
ten  in  the  following  form: 


to  nodal  locations.  Using  the 
inverse  of  Equation  2.1  is  writ- 


(  \ 

r 

Cl 

1 

2A 

ai 

a2 

a3 

L1 

■  C2 

bi 

b2 

b3 

- 

L2 

c3 

V.  J 

ci 

C2 

C3 

.  L3- 

2.3 


Substituting  Equation  2.3  into  2.1  yields: 
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L  =  2A  [1V2] 


Letting: 


ai  a2  a3 
bl  b2  b3 


C1  C2 


3 

%  / 


2.4 


[a]  = 

[ala2a3] 

[b]  = 

[bib2b3] 

[c]  = 

[C1C2°3] 

2.5 


<Li> 


L3J 


Substituting  2.5  into  2.4  results  in: 

L  ■  2X  Cta]  +  tblxi  +  [c]x2){L.} 


2.6 


Referring  to  Equation  1.1  it  is  apparent  that  when  is  a  constant, 


u„ 


'3  =  Si  +  ?2X1  +  ?3X2 

Following  a  similar  procedure  to  that  above  results  in: 


2,7 


u„ 


2X  +  [b]x  +  [c]x2){u3  } 

i 


2.8 


where  {u^.}  are  nodal  displacements  of  node  i  in  the  X3  direction.  The 
virtual  work  of  the  distributed  load  is: 


dW¥ 


du^  L  d  A 


2.9 


Substituting  Equations  2.6  and  2.7  into  2.9  results  in: 


dW  = 

Lt 


A  d  {u3  +  [bl Xl+ [c] X2)T2i( [a] + [b] xi+ [£] x2^ {Li }dA  2,10 
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The  virtual  work  done  by  the  distributed  load  L  is  now  replaced  by  effec¬ 
tive  nodal  forces  which  do  virtual  work  as  given  by: 

dW  =  d  {u3  }  {p  }  2.11 

i  E 

where  Wp  represents  the  work  done  by  the  effective  concentrated  forces 
Fg  .  Since  the  effective  forces  are  caused  by  the  distributed  loads, 
it  follows  that: 

dW_  =  dWT  2.12 

E  L 

Substituting  Equations  2.10  and  2.11  into  2.12  gives: 


{V 


4?  w  {Li> 


2,13 


where : 


[D]  =  [a]T  [Dj]  +  [b]T  [D2]  +  [cf  [D3] 


[Dj]  =  [a]  +  [b] 


AxjdA  +  [c] 


AX2dA 


[D21  =  [a] 


AxxdA  +  [b] 


AX2  dA  +  [C] 


AXlX2dA 


2.14 


[D3]  =  [a] 


Ax2dA  *  [b] 


AX!X2dA  +  [c] 


Ax2dA 


The  reason  for  this  method  of  writing  the  final  result  is  that  this 
method  was  found  computationally  more  efficient  to  calculate  the  final 
forces.  Section  6  will  discuss  the  loading  used  in  the  test  examples. 
Knowing  how  to  calculate  the  internal  and  external  forces,  the  next 
step  is  the  solution  of  the  nodal  equilibrium  equations. 


23 


III.  DYNAMIC  ANALYSIS 


Because  the  use  of  a  variable  time  increment  would  be  advantageous 
in  analyzing  a  structure  with  a  load  applied  over  a  specific  length  of 
time,  a  method  of  analysis  was  developed  to  include  the  variable  time 
increment.  This  gives  the  programmer  the  responsibility  to  choose  appro¬ 
priate  time  increments  over  the  length  of  time  that  the  structure  is 
analyzed.  Good  engineering  judgment  will  dictate  the  size  of  these 
time  increments.  When  there  is  a  rapid  change  in  external  load,  the 
time  increments  should  be  small  when  compared  to  times  when  there  is  less 
change  in  load  or  no  load  at  all.  Normally  the  initial  analysis  should 
be  very  accurate  so  even  smaller  time  increments  could  be  used.  The 
numerical  solution  of  Equation  1.59  by  the  finite  difference  method  in¬ 
volves  the  replacement  of  the  time  derivatives  by  their  finite  differ¬ 
ence  equivalents . 2 > 3  p0r  convenience  the  matrix  notation  will  be  drop¬ 
ped  and  the  displacement  symbol  (A)  will  be  replaced  by  the  quantity  x 
with  a  subscript  defining  the  time  interval.  In  this  approach  the  time 
history  is  divided  into  discrete  time  intervals  whose  length  will  be  de¬ 
noted  by  h  with  a  subscript  defining  the  time  interval.  From  the  gen¬ 
eral  theory  of  kinematics,  the  velocity  and  displacement  relations  are 
written  for  the  n^1  and  (n+l)^  time  interval  as  follows: 


x 

n 


x  +  ^-  (x  +  x  ) 
n-1  2  n-I  n 


3.1 


x 

n 


n-1 


+  h  x  T  +  (4  -  B)h2x  +  Bh2x 
n  n-1  2  n  n-1  n  n 


3.2 


n+1 


n 


+  h 


n+I  n 


+  (t  -  3)h 


n+1  n 


+  Bh 


n+1  n+1 


3.3 


where  x,  x  and  x  represent  acceleration,  velocity,  and  displacement  at 
time  intervals  denoted  by  the  subscript. 

Equation  3.1  means  that  the  velocity  at  the  end  of  the  interval  is 
equal  to  the  sum  of  the  velocity  at  the  beginning  of  the  interval  and 
the  product  of  the  time  interval  and  the  average  of  the  acceleration  at 
the  beginning  and  end  of  the  interval.  Equation  3.2  and  similarly 
Equation  3.3  are  obtained  by  integrating  Equation  3.1  and  introducing  a 
weighted  acceleration  parameter,  3,  to  express  the  average  acceleration. 
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For  example,  3=1/4  gives  a  linear  estimation  for  the  average  accelera¬ 
tion.  The  coefficient  3  is  chosen  so  as  to  best  represent  the  system 
being  analyzed.  Newmark^  discusses  convergence  and  the  values  of  beta 
and  shows  that  for  dynamic  problems,  3=1/4  gives  an  infinite  stability 
limit  and  is  the  value  used  for  most  dynamic  problems.  At 

n+1  n  n-1 

t  =  l  h,  ,  £  h,  ,  £  h,  respectively,  the  equations  of  motion  become: 
k=l  k  k=l  k  k=l  k 


M  x 


n+1 


n+1 


3.4 


M  x  =  F 
n  n 


3.5 


M  x 


n-1 


n-1 


3.6 


F 


FI  +  fe 


3.7 


The  unknowns  in  Equations  3.1  through  3.6  are  xn+j,  x  ,  xn-i,  xn+i,  xn, 
xn-i-  Equations  3.3  will  be  used  to  solve  for  displacement  at  time 
n+1.  Equations  3.4,  3.5,  3.6  give  values  of  the  accelerations  in  terms 
of  known  forces  and  masses.  The  only  unknown  then  is  xn.  Rearranging 
Equation  3.1  gives: 


.  un  11  u  n  n 

x  =  x  .+•— -  x  +7^-  X  +3h  x  -3-^-  X  .  +  3;^-  X  -3h  x  3.8 

n  h  n-1  2h  n-12  n  n  n-1  Hi  n-1  h  n  nn 
n  n  n  n 


This  can  be  rewritten  as: 


x  =  [h  x  +  -8)h2x  +8h2x  ]  +  8h  x  T  +  (^-  -3)h  x 

n  h^  L  n  n-1  ^2  n  n-1  n  nJ  n  n-1  ^2  J  n  n 


3.9 


Combining  Equations  3.2  and  3.9  results  in  the  following: 


*  1  . . 

x  =  7—  (x  -x  T)  +  8h  x  .  +  f—  -  3)h  x 
n  h  v  n  n-ly  K  n  n-1  ^2  J  n  n 

n 


3.10 


Inserting  Equation  3.10  into  3.3  yields: 


25 


3.11 


x  =x  +-r — (x  -x  t)  +  (h  t+h  )h  , (—  -3)x  +  3h2  ,x  +3h  ,h  x  T 
n+1  n  h  n  n-r  ^  n+1  n+1  ^2  J  n  n+1  n+1  n+1  n  n-1 

n 


Substituting  Equations  3.4,  3.5,  3.6,  3.7  into  3.11  gives: 


x  =  x  4 — -(x  -x  )-  (h  _+h  )h  T(^  -m  1(Ft  +Fc  ) 
n+1  n  h  v  n  n-r  v  n+1  n  n+1 ^2  y  1  E 

n  n  n 


,2  „-l~  ,  w-1 


'F1  /FE  ,)+Sh„hn*lM  tpl  +FE  ’ 
n+1  n+1  n-1  n-1 


3.12 


Using  the  original  matrix  notation  this  becomes: 

{A}  =  {A}  +-£ii({A}  -{A}  T)+h  _(h  .+h  ) 4  -3) [M]_1({F.}  +  {F  }  ) 

n+1  n  h  v  n  n-1  n+lv  n+1  nJ  K2  1  n  E  n 

n 


*eh^1tM]-1({FI)nilHFE}n,1).ehntlhntMrI({F1)n.1.{FE}n.1) 


3.13 


For  this  method  to  work,  it  is  assumed  that  the  internal  forces  are 
known.  But  to  know  the  internal  forces,  the  stresses  must  be  known  as 
indicated  in  Equation  1.58.  To  know  incremental  stress,  the  strain  must 
be  known.  To  know  strains,  deflections  must  be  known.  If  the  time 
steps  are  small  enough,  then  there  is  not  much  change  in  stress  from  one 
time  increment  to  another.  Thus,  this  method  is  extremely  accurate. 

But  for  larger  time  increments  another  approach  must  be  used.  This  is 
described  in  Sections  4  and  6. 
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IV.  MODIFICATION  FOR  LARGE  TIME  INCREMENTS 


The  limitations  on  the  size  of  the  integration  time  step  are  a  di¬ 
rect  result  of  the  finite  element  model  which  has  a  lumped  mass  at  each 
node.  As  the  numerical  integration  proceeds,  the  masses  move  relative 
to  each  other.  If  the  time  step  is  too  large,  the  relative  motion  of 
the  masses  is  exaggerated  and  artificial  oscillation  is  induced.  In 
order  to  combat  this  artificial  oscillation,  the  x^-X2  deflections  are 
coupled  and  the  deflections  in  the  X3  direction  are  coupled  through  the 
thickness. 

First,  then,  the  in-plane  (u^-^)  displacements  will  be  discussed. 
In  order  to  have  a  more  convenient  equation  to  work  with,  Equation  3.13 
is  written  with  a  constant  time  increment  as  follows: 

{A)n+1  =  2{A}  -{A}  +0h2[M]'1[{F  }  +d-2){F  }  +{F  }  ] 

1  n+1  p  1  n  n-1 

+£h2 [M] _1 [{F  }  +(I  _2){Fp}  +{Fp}  ]  4.1 


where  {A}  is  the  displacement  matrix,  3  is  the  acceleration  parameter, 
h  is  the  time  interval,  {Fj }  is  the  internal  force  matrix,  {Fg}  is  the 
external  force  matrix,  [M]  is  the  mass  matrix,  and  the  subscripts  n, 
n-1,  and  n+1  denote  time  intervals. 

In  analyzing  the  and  U2  displacements,  it  is  assumed  that  they 
are  linearly  dependent  through  the  thickness.  This  forces  plane  sect¬ 
ions  to  remain  plane. 

This  first  assumption  results  in  the  following  equations: 


U1  =  qi  +  zq2 


U2  =  q3  +  Zq4 


4.2 


where  q^,  k=l,  4,  are  unknown  generalized  coefficients  called  the  trans¬ 
formed  displacements  and  z  is  the  distance  in  the  x^  direction  of  the 
node  from  the  center  of  gravity.  The  importance  of  having  z  be  the  dis¬ 
tance  from  the  center  of  gravity  will  be  discussed  when  the  transformed 
mass  matrix  is  discussed.  In  matrix  notation  Equation  4.2  becomes: 
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{A}  =  [tf] {q} 


4.3 


where  {q}  is  the  matrix  of  transformed  displacements  and  [tf]  is  the 
transformation  matrix  described  below.  Letting  Z  be  the  number  of 
layers  of  material  and  i=£,  l+£  be  the  nodal  location  in  the  thickness 
direction,  the  transformation  matrix  can  be  written  as: 


4.4 


where: 


1  z.  0  0 


i 


4.5 


0  0  1  z. 


l 


and  z^  is  the  distance  of  node  i  from  the  center  of  gravity.  The  nota¬ 
tion  i=l  is  the  node  at  the  bottom  of  the  plate  and  i=l+£  is  the  node  at 
the  top. 

Since  the  displacements  are  written  in  terms  of  transformed  dis¬ 
placements,  the  forces  should  be  written  in  terms  of  the  transformed 
forces.  Letting  {fg}  be  the  matrix  of  external  forces  corresponding  to 
{q},  and  {fj}  be  the  internal  forces  also  corresponding  to  {q},  the 
principle  of  virtual  work  states : 


d  {A}T  {Fj}  =  d  {q}T  {f^ 


4.6 


Transposing  Equation  4.3  yields: 


4.7 


Substituting  Equation  4.7  into  4.6  yields: 


d  {q  }T  [tf  ]T  {Fj }  =  d  {q  }T  {fj} 


4.8 
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Therefore; 


{fj}  =  [tf]T  {Fj} 

4.9 

Similarly: 

{fEl  =  [tf]T  {fe} 

4.10 

Finding  a  transformed  mass  matrix 

[m]  starts  with: 

-[M]  {A}  =  {F}. 

L  J  inertia 

4.11 

By  virtual  work: 

-d  {A}T  [M]  {A}  = 

d{A}T  {F}.  ,. 

inertia 

4.12 

-d  {q}T  [m]  {q}  = 

d^>T  <F>inertia 

4.13 

Therefore : 

d  {A}T  [M]  {A}  = 

d{q}T  [m]  {q} 

4.14 

Substituting  Equations  4.3  and  4.7 

into  4.14  yields: 

d  {q}T  [m]  {q}  = 

d{q}T  [tf] T [M]  [tf]  {q} 

4.15 

Dividing  out  the  unnecessary  terms 

gives: 

[m]  =  [tf]T  [M]  [tf] 

4.16 

Now  the  reason  for  being  the  distance  from  the  center  of  gravity 
will  become  apparent.  The  calculations  are  much  simplified  by  the  mass 
matrix  being  a  diagonal  matrix  as  in  the  case  for  the  original  mass  ma¬ 
trix.  The  original  mass  matrix  was: 

[M]  =  [M  ]  4.17 

Performing  the  matrix  multiplication  in  Equation  4.16  using  Equations 
4.4,  4.5  and  4 . 7  produces : 


lt£ 

l  M. 

i=£  x 

l+£ 

y  M.  z. 
i=£  1  1 

0 

0 

l+£ 

y  M.z. 
i=£  1  1 

l+£ 

y  M.z? 

u  IT 

i=£ 

l+£ 

y  m.z. 

i=£  1  1 

0 

II 

£ 

0 

l+£ 
y  M.z. 
i=£  1  1 

l+£ 

1  Mi 

i=£ 

l  +  £ 

y  m. 

i=£  ■ 

0 

0 

l+£ 
y  M.z. 
i=£  1  1 

l+£ 

y  m. 

i=£ 

But, 

l+£ 

y  M.z.  =  0  4.19 

i=£  1  1 

by  the  definition  of  the  center  of  gravity,  thus  causing  [m]  to  be  a 
diagonal  matrix. 

Referring  to  Equation  4.1,  the  only  elements  that  still  must  be 
transformed  are  the  deflections  for  past  time  intervals.  From  Equation 


{q}  =  [tf]"1  {A}  4.20 
n  n 

To  find  these  generalized  displacements,  it  is  only  necessary  to  know 
four  of  the  actual  displacements  to  solve  for  four  unknowns.  This  is 
easily  done  by  hand  and  results  in: 


[tf:] 

[tf2] 


VZ2 


zrz2 


zrz2 


zrz2 


zrz2 


zrz2 


zrz2 


zrz2 


4.21 
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Substituting  the  transformed  quantities  into  Equation  4.1  yields: 


{q}n+i  =  2{q}n-{q}n_1+Bh2[m]’1 


{fj}  + 

b2 

<V  /{fi}  . 

n 

b  j 

n-1  n-2 

Ww1 


n+l 


'i-2' 


n  n-1 


4.22 


It  should  be  noted  that  the  internal  transformed  forces  are  displaced 
by  one  time  increment.  Because  these  forces  are  small,  even  at  large 
time  increments,  this  yields  accurate  results  with  small  error.  Equa¬ 
tion  4.22  is  thus  solved  and  Equation  4.3  transforms  the  results  to 
global  displacements. 

Although  this  time  lag  is  acceptable  for  the  in-plane  displacements, 
it  is  not  acceptable  for  the  U3  displacements.  The  reason  for  this  is 
that  the  external  force  is  being  applied  in  the  U3  direction,  thus 
making  these  internal  forces  more  sensitive  to  larger  time  intervals. 

In  order  to  account  for  the  change  in  the  internal  force,  a  model  was 
sought  to  couple  the  deflections  through  the  thickness. 

In  finding  a  model  to  represent  what  happens  through  the  thickness, 
it  is  necessary  to  see  what  the  unknowns  are.  From  Equation  4.1  the  un¬ 
knowns  are  {A)n+i  and  {F j }n+i -  All  the  other  terms  are  known.  In  order 
to  predict  what  {Fi}n+i  is,  it  is  necessary  to  couple  the  deflections 
through  the  thickness  and  to  assume  all  strains  small  when  compared  to 
the  strain  in  the  u^  direction.  This  can  be  done  by  letting: 


{F  }  =  {F  }  +  {AF  } 

n+l  n  n+l 


4.23 


where  {AFj}n+j  is  the  change  of  the  internal  force  between  time  inter¬ 
vals.  The  deflection  in  the  U3  direction  is  then  coupled  by  the  model 
shown  in  Figure  5.  This  model  assumes  the  stiffness  between  the  nodal 
points  in  the  thickness  direction  is  much  greater  than  the  stiffness  be¬ 
tween  in-plane  nodal  points.  As  long  as  the  external  force  is  in  the 
u^  direction  this  is  a  good  assumption. 

From  Figure  5: 


AFt  =  k.[(A.  .  .-A.  J-CA.  .  -A.  )]  -  k.  . 

I.  ,  1  l+l, n+l  i,n+l  l+l, n  i,n  l-l 

i,n+l 

[(A.  ,  -A.  .  J-(A.  -A.  .  )] 

LL  i,n+l  l-l, n+l  i,n  i-l,n  J 


4.24 
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where  i  refers  to  the  nodal  location  through  the  thickness  and  n  refers 
to  the  time  increment.  The  predicted  stiffness  (k^)  is  found  from  the 
orthotropic  properties  (C^j ,  i=l,6,  j=l,6).  In  matrix  notation: 


f  \ 

all 

"Cll 

C12 

C13 

0 

0 

0 

t  y 

en 

°22 

C21 

C22 

C23 

0 

0 

0 

e22 

*°33 

,  = 

C31 

C32 

C33 

0 

0 

0 

- 

e33  ' 

ai2 

0 

0 

0 

C44 

0 

0 

E12 

a23 

0 

0 

0 

0 

C55 

0 

e23 

°13. 

0 

0 

0 

0 

0 

C66 

.e13- 

Using  the  assumption  that  all  strains  are  small  when  compared  to  the 
strain  in  the  u^  direction  yields: 

a33  =  C33e33  4.26 


This  gives  the  stiffness  for  a  unit  cube  equal  to  C33. 


k. 

1 


C33.Areai 

1 


t. 

1 


Therefore: 


4.27 


where  t^  is  the  thickness  of  layer  i,  Area^  is  the  area  used  to  compute 
the  mass  of  the  nodal  point,  and  C33.  is  the  orthotropic  property  of  the 

quadrilateral  that  the  node  lies  in.  Letting: 


A7  =2A.  +A.  + 

i,n+l  i,n  i,n-l 


(3h 


M. 

1 


1  0 

+ 

-2 

I. 

3 

i,n 

j 

FI.  +FI.  , 
i,n  i,n-l 


M. 

1 


Ei,n+1 


r  „ 1  > 

+ 

1-2 

F_  +F_ 

8 

E .  E .  . 

J 

i,n  i,n-l_ 

and  substituting  Equations  4.23  and  4.28  into  4.1  produces: 
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Ai,n+1  =  Ai,n+1  +  ^M.  *-ki  ^i+l,n+l"Ai,n+l^  "  ki-l  ^i,n+l_Ai-l,n+l^ 


k.  (A.  .  -A.  )  +  k.  .  (A.  -A.  ,  )] 

1  l+l, n  i,n  l-l  i,n  i-l,n 


4.29 


where  the  only  unknowns  are  the  deflections  at  time  n+1.  This  produces 
i=l+£  (£  is  the  number  of  layers)  number  of  simultaneous  equations 
which  can  be  solved  for.  The  technique  of  coupling  the  motion  of  the 
masses  in  the  thickness  direction,  as  described  in  this  section,  per¬ 
mits  the  use  of  larger  time  increments  than  otherwise  possible.  The 
effect  of  this  coupling  is  to  reduce  the  numerical  oscillations  produced 
by  excessive,  relative  motion  of  the  adjacent  masses.  However,  the  use 
of  larger  time  increments  introduces  an  additional  side  effect  for 
rapidly  varying  external  loads.  If  these  loads  vary  by  a  great  amount 
from  one  time  interval  to  another,  then  the  transverse  stresses  from  the 
n^  time  interval  will  be  greatly  out  of  balance  with  the  external 
forces  at  (n+l)^*1  interval.  This  leads  to  large  force  unbalance  which 
can  lead  to  numerical  errors.  The  Section  VI  of  the  report  will  dis¬ 
cuss  a  method  for  eliminating  these  large  differences. 
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V.  PLASTICITY 


Orthotropic  Elastic-Plastic  Yielding 

The  present  analysis  allows  for  permanent  deformation  to  occur  in 
the  structure.  Because  of  the  practical  considerations  each  finite  ele¬ 
ment  is  assumed  to  be  either  elastic  or  plastic.  The  general  method  of 
attack  is  to  calculate  the  stresses  as  if  they  were  elastic.  Using 
Hill's  orthotropic  yield  criterion, ^  the  stresses  are  checked  to  see  if 
they  are  physically  compatible  with  the  yield  criterion.  If  they  are 
compatible  the  analysis  continues  with  those  values  of  stress.  If  they 
are  not  compatible,  then  plastic  flow  has  occurred  and  using  the  flow 
rule  and  the  yield  criterion,  the  stresses  are  recalculated  to  account 
for  this  plastic  behavior. 

The  first  step  then  is  to  calculate  the  stresses.  Using  tensor 
notation,  the  incremental  form  of  the  stress-strain  equation  for  elas¬ 
tic  case  gives  the  stress  change  from  time  tn  to  tn+^: 


C.  .,  nde, 
ljkl  kl 


5.1 


where  the  superscript  T  indicates  a  test  value.  The  total  stress  at  a 
time  (tn+1) : 


aT. 

ij 

where  the  stress  a 


=  a.  . 


n 


T 

da  .  . 
il 


5.2 


il 


is  from  time  t 


n 


n 


This  value  is  then  put  into  Hill's  yield  criterion  which  requires  six 
yield  stresses  (Yjj).  Before  showing  the  yield  criterion,  the  following 


constants  are  defined: 


11 


11 


22 


33 


22 


22 


11 


33 


34 


5.3 


33  2  2  2 

33  11  22 


Then  the  yield  criterion  can  be  written  as: 


2  2  2  2  2  2 

offer  i-^  +  f22_  °33  °12^23  ^13  -  - 

aij  y2  y2  2  2  2  *2  +  Yll°22°33  +  Y22°lia33 

11  22  33  12  23  13 


+  Y33aHa22  1 


5.4 


T  T 

Using  Equation  5.2  in  5.4,  and  if  2f(a. . )  <^  1,  then  a. .  represents 


il' 


il 


the  real  stress  at  t  but  if  2ffO. .)  >  1,  yield  has  occurred. 

n+1  l  y 

To  calculate  the  stress  for  plastic  flow,  the  strain  increment  is 

6  P 

divided  into  elastic  (e  )  and  plastic  (e  )  strain. 


de.  . 
il 


de?.  +  de?. 
il  il 


5.5 


The  flow  rule  is  written  in  the  following  manner: 

On 


de^  =  d\ 


'11 


Y2 

l.  11 


Y22a33+Y33°22 


=  dAT 


11 


de22  "  d* 


'22 


22 


Ylla33+Y33CT11^ 


dAT 


22 


de33  =  dX 


33 


33 


Y  a  +Y  a 
11  22  22  11 


=  dAT 


33 


de 


P  _ 
12 


dA 


12 


12 


=  dAT 


12 
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5.6 


dE23  =  dX 


23 


23 


=  dAT 


23 


de?„  =  dA  -P-  =  dAT 


'13 


13 


13 


where  Tji  are  defined  by  Equations  5.6  and,  for  the  purpose  of  numeri¬ 
cal  calculations,  they  are  assumed  to  be  given  by  the  stresses  from 

time  t  .  In  short: 
n 


de?.  =  dAT.. 
iJ  iJ 


5.7 


In  view  of  Equation  5.5,  the  actual  form  of  Equation  5.1  is: 


da..  =  C.  ...def, 
ij  ijkl  kl 

Substituting  Equation  5.5  into  5.6  results  in: 

doij  ’  WdEki  -  dEU> 


5.8 


5.9 


Now  Equations  5.1  and  5.7  are  substituted  into  5.9  giving: 


da. .  =  da.  .  -  C.  ..  .dAT, . 

ij  ij  ljkl  kl 


5.10 


Using  Equation  5.10  to  update  the  stress  at  t  gives  stress  at  t^+ ^ : 


aij  aij  ‘  dXCijklTkl 


5.11 


Letting  T^.=  Equation  5.11  becomes: 


aij  aij  dXTij 


5.12 


Substituting  Equation  5.12  into  5.4  results  in  a  quadradic  equation  in 
dx  as  follows : 


Ad  A  -  2Bd  A  +  C  =  0 


5.13 


where : 
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YTT+YTT+YTT 
11  22  33  22  11  33  33  22  11 


B  = 


TliaIl  T22°22  ,  T33°33  T12<4  T23°23  T13°l3 


11 


22 


33 


12 


23 


13 


f  \ 

—  T  —  T 

V  T  -  T 

f  \ 

—  T  -  T 

T22a33+T33a22 

+  Y 

Tll033+T33ail 

+  Y 

Tlia22+T22ail 

2 

V  J 

22 

2 

^  J 

33 

2  J 

C  =  2f faT.)  -  1 

ij 


5.14 


Solving  Equation  5.13  results  in: 


,,  _  2B±/(2BT^4AC 

d  2A  5.15 

Because  C  is  defined  in  Equation  5.14  as  the  yield  criterion,  as  C 
approaches  0,  dA  approaches  0.  Clearly,  the  minus  sign  is  the  only  sign 
that  is  physically  acceptable.  Multiplying  top  and  bottom  of  Equation 
5.15  by  B+/B*-AC  produces: 


dA 


C 

B+/B^-AC 


5.16 


This  dA  is  substituted  into  Equation  5.12  to  give  the  actual  stress 
at  time  tn+i. 

Orthotropic  Elastic-Viscoplastic  Yielding 

Because  of  the  nature  of  an  impulsive  load  acting  on  a  material,  the 
plastic  yielding  is  such  that  it  is  history  dependent.  To  be  able  to 
handle  materials  with  viscous  coefficients,  a  viscoplastic  analysis  has 
been  developed. 
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In  this  analysis  it  is  assumed  that  the  strain  is  divided  into 
elastic  (dee)  and  viscoplastic  (devP)  strain. 


de.  .  =  de? .  +  deY? 

Using  Hill's  Flow  Rule,  Equation  5.7  is  modified  to  be: 


deY?  =  dXT.  . 


13 


5.17 


5.18 


It  should  be  noted  that  the  viscoplastic  strain  changes  satisfy  the  in¬ 
compressibility  condition: 


l  deY? 

i=l 


11 


5.19 


The  quantities  T^j  represent  six  independent  quantities  and  they  can  be 
arranged  in  a  matrix  form  and  then  can  be  related  to  a  stress  matrix  as 
follows : 


{T}  =  [H] {a} 

where : 


0 


0 


0 


0  0 

0  0 

0  0 


0  0  0 

0  0  0 

0  0  0 


10  0 


5.20 


5.21 
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By  comparing  the  flow  rule  in  Equation  5.18  to  an  isotropic  case,  it  can 
be  noted  that  the  quantities  T^j  have  the  same  role  as  the  deviatoric 
stresses  and  e^P  strains  as  the  deviatoric  strains.  In  fact,  it  may  be 
noted  that  Equation  5.18  reduces  in  the  limit  to  the  isotropic  case  when 
the  material  properties  are  the  same  in  all  three  directions.  Conse¬ 
quently,  following  the  procedure  developed  for  the  Bingham  material, ^ 
the  strain  rate  dependence  is  introduced  by  defining  {TF}: 


{TF}  =  {T}  +  n{eVp}  5.22 

where  q  represents  viscous  coefficient,  and  {T}  is  the  quantity  which 
satisfies  the  yield  criterion.  By  using  Equation  5.20,  a  final  stress  is 
defined  as: 


{aF}  =  [H]  1  {TF}  5.23 

If  no  yielding  has  occurred,  the  viscoplastic  strain  increment  is  zero. 
So  this  analysis  begins  with  a  trial  incremental  stress  where: 

{daT}  =  [C]  {de}  5.24 

and  [C]  is  the  orthotropic  relationship  between  stress  and  strain. 
Equation  5.24  is  inserted  into  Equation  5.2  and  this  value  is  inserted 
into  the  yield  criterion  in  Equation  5.4.  If  yield  does  not  occur,  the 
trial  stress  is  equal  to  o^.  However,  if  yieldpdoes  occur  then  Equa¬ 
tions  5.22  and  5.23  must  be  used  to  calculate  o  as  follows: 

From  Equation  5.17  and  5.24: 

{daF}  =  [C]  {de}- [C]  {deVp}  5.25 

Then  as  in  Equation  5.2: 

{aF}  =  {aT}  -  [C]  {de^}  5.26 

Multiplying  Equation  5.22  by  [H]"1  and  using  Equation  5.23  and 
5.24,  one  gets: 


{/}  =  {a}  +  n  [H] -1  {e^}  5.27 

From  Equation  5.26  and  5.27,  when  solving  for  {a}  one  finds: 


{a}  =  {aT}  -  [C]{devp}-n[H]"1{evp} 


5.28 
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By  using  Equation  5.28  it  is  easily  shown  that; 

{eVp}  =  Wevp}  =  ~  {T}  5.29 

Substituting  5.18  and  5.29  into  5.28  produces: 

{a}  =  {cT}  -  dX[C]{T}  -  5.30 

Defining  another  variable  {T } : 

{T}  =  [C]  {T}  5.31 

and  using  the  inverse  of  Equation  5.20  and  Equations  5.31  in  5.30  yields 


{a} 


{a  }  -  dX 


m  *  {a  j 


5.32 


The  question  in  a  dynamic  problem  always  arises  as  to  what  value  of 
stress  is  used  for  the  flow  rule.  In  this  formulation  the  stress  used 
in  the  flow  rule  is  approximated  by  the  trial  stress.  A  closer  approxi¬ 
mation  can  be  formed  by  doing  an  iterative  loop  on  this  equation,  but 
little  difference  is  found  in  the  solution  when  this  is  done. 

The  {o}  stress  formed  here  in  Equation  5.32  is  input  in  the  yield 
criterion  and  results  in  the  following  equation: 


Ad  A  -  2BdX  +  C 


5.33 


where : 


t:  . 


T1  +  —  aT 
ij  At 


5.34 


rp  •  2  rp  •  2  rp  *  2  rp  •  2  rri*2 

,  11,  22,  33,  12  23  13  y  t'T"  +  t  *  t  "  +  T  ’  T  " 

a  =  t 'zrzrzrzrzr  +YnT22T33+Y22TnT33+Y33T22Tn 


ii 


22  33 


12 


23  13 


5.35 


rp  rp  rp  rp  rp  rp 

T"  rr1  t"  a  T‘  a1  T‘  a1  T‘  a1  T'  cr1 
11  11  22 22  1 33°33  X12  12  23  23  13  13 

B  =  -  +  -  +  -  +  -  +  -  +  - 


11 


Y2 

22 


Y2 

^>3 


12 


23 


13 


40 


'  •> 

T22033+T33a22 

4-  Y 

Til°33*T33°n 

4-  Y 

T  T 

T"  a  +T*  a 

11  22  22  11 

2 

V.  J 

Y22 

2 

k.  J 

Y33 

2 

5.36 


C  =  2  f  (aT.)  -  1 

1J 


5.37 


where  2f(a?^.)  =  l  is  the  yield  function. 

The  proportionality  constant  dX  can  be  solved  for  and,  as  shown  for 
Equation  5.13,  is: 


B+/BZ-AC  5.38 

This  value  of  dX  is  then  substituted  into  Equation  5.32.  By  using 
Equations  5.28,  5.20  and  5.27: 


{/}  = 


1  + 


qdX 


At 


{a} 


5.39 


With  this  analysis  completed,  the  problem  of  the  oscillating 
stress  is  reduced  but  still  sometimes  occurs.  This  then  leads  to  the 
next  section  which  deals  with  the  element  equilibrium  equations. 


41 


VI.  ELEMENT  EQUILIBRIUM  EQUATIONS 


The  previous  analysis  helped  allow  for  the  use  of  larger  time 
increments  in  the  numerical  time  integration.  Although  the  results 
showed  the  deflection  to  be  stable,  there  were  normal  stress  oscilla¬ 
tions.  In  order  to  compensate  for  these  small  errors  in  strain,  the 
external  pressure  is  used  to  develop  element  equilibrium  equations. 

The  modification  of  the  analysis  involves  the  introduction  of  addi¬ 
tional  dynamic  equations  which  express  the  transverse  acceleration  of 
the  elements  as  a  whole,  in  addition  to  the  nodal  accelerations  as  ex¬ 
pressed  by  Equation  1.56.  Consider  a  small  quadrilateral  section  of 
the  plate  as  illustrated  in  Figure  6.  In  general,  this  section  will  be 
composed  of  N  layers  and  loaded  by  external  distributed  force 
p(xj,X2,t).  The  transverse  acceleration  of  this  plate  section  that  is 
accelerating  in  the  X3  direction,  is  produced  by  two  distinct  forces. 

The  first  of  these  is  the  external  force  p(xi,X2,t),  and  the  second  is 
the  bending  force  in  the  plate.  As  a  matter  of  observation,  it  is  logi¬ 
cal  to  expect  the  external  force  to  be  predominant  in  the  early  stages 
of  the  time  history  while  the  bending  forces  should  become  more  impor¬ 
tant  as  the  deformation  of  the  plate  increases.  Since  the  amount  of 
each  contribution  will  vary  with  time,  the  relative  amounts  of  acceler¬ 
ation  from  each  load  will  become  one  of  the  unknowns  of  the  problem. 

In  order  to  take  this  into  consideration,  it  is  assumed  that  the  amount 
of  acceleration  due  to  the  external  force  will  be  related  to  the  total 
average  acceleration  by  a  relation: 

a  =k  a  6.1 

P  t 

where  ap  is  the  acceleration  due  to  the  force  p(xj,X2,t),  at  is  the 
total  acceleration,  and  k  is  an  unknown  parameter.  It  may  be  noted 
that  the  accelerations  in  Equation  6.1  will  vary  through  the  plate  thick¬ 
ness.  Consider  now  the  acceleration  in  the  X3  direction  of  the  indivi¬ 
dual  layers  in  the  plate  sections  of  Figure  6.  These  layers  are  now 
the  actual  finite  elements.  The  normal  stress  in  the  X3  direction  will 
vary  from  layer  to  layer,  and  the  top  and  bottom  values  of  this  stress 
will  be  equal  to  the  applied  surface  loads.  The  plate  elements  through 
the  thickness  of  the  plate  are  illustrated  in  Figure  7.  The  interface 
stresses  are  denoted  by  a  with  an  appropriate  subscript.  Writing  the 
dynamic  equations  for  each  layer  results  in  the  following  N  equations: 
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a2-°I  =  V1'1 


°3  -  °2  =  ap2P2t2 


6.2 


°N  ”  °N-1 


.  -  a,T  =  a  p,Tt,T 
N+l  N  pXT  N  N 

N 

where  N  represents  the  total  number  of  layers,  is  the  material  density 
and  t  is  the  thickness  of  each  layer.  It  may  be  noted  from  the  boundary 
values  that: 


a  =  0  6.3 

Vi  =  p^i-v5 


Using  Equations  6.1  and  6.3  and  adding  Equation  6.2  results  in  the 
solution  of  the  unknown  constant  k: 


k 


p(x1,x2,t) 


N 


l 

3=1 


a  p.t. 
t.  j  i 
3 


6.4 


Since  the  finite  element  is  formulated  in  terms  of  element  stresses 
rather  than  the  interface  stresses,  the  latter  have  to  be  expressed  in 
terms  of  the  element  stresses.  This  is  done  by  linear  extrapolation: 


a. 

3 


(a.+a.  . ) 

3  3-1 


6.5 


where  aj  represents  the  stress  033  (normal  stress  in  the  X3  direction) . 
By  using  Equations  6.4  and  6.5  in  Equation  6.2  it  is  possible  to  solve 
for  the  element  stresses  to  cN.  This  procedure  will  yield  the  normal 
stress  in  the  X3  direction  for  each  element  in  the  plate. 
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In  applying  this  analysis  to  the  existing  finite  element  model,  it 
must  be  realized  that  Equations  6.2  are  not  completely  new  dynamic  equa¬ 
tions;  but  rather  the  information  contained  in  these  equations  is  al¬ 
ready  contained,  in  a  different  form,  in  the  original  dynamic  equation 
represented  by  Equation  1.59.  However,  Equations  1.59  represent  nodal 
accelerations,  but  the  accelerations  in  Equations  6.2  are  some  average 
accelerations  in  the  elements  and  therefore  the  two  sets  of  accelera¬ 
tions  are  related  but  are  different  quantities.  In  fact,  in  the  present 
analysis  the  element  accelerations  are  defined  by  averaging  the  nodal 
values.  In  order  to  see  how  Equations  6.2  are  incorporated  into  the 
analysis,  it  is  useful  to  review  how  the  original  analysis  proceeded. 
Equations  1.59  are  solved  by  a  finite  difference  procedure  and  the  re¬ 
sulting  displacements  are  used  to  evaluate  strains  and  therefore 
stresses.  The  stresses  are  then  used  to  check  for  yield  and  to  evaluate 
the  body  forces  for  the  next  time  interval .  In  the  modified  approach 
Equations  1.59  still  solve  for  the  displacements  which  are  used  to 
evaluate  all  the  strains.  The  displacements  from  Equations  1.59  are 
also  used  to  evaluate  the  element  accelerations  which  are  needed  for 
solution  of  Equations  6.2.  The  element  accelerations  are  related  to 
the  nodal  displacements  by  averaging  over  all  the  corner  nodes  for  each 
el ement . 


The  value  of  the  normal  stress  as  calculated  by  this  method  is 
stored,  and  the  previous  values  of  stress  are  used  to  solve  for  the 
plastic  yielding.  After  the  completion  of  the  plastic  yielding,  the 
stresses  are  modified  so  that  they  remain  on  the  yield  surface,  yet  com¬ 
ply  with  the  element  equilibrium  equations.  This  is  easily  accomplish¬ 
ed  by  referring  to  Hill's  yield  criterion  (Equation  5.4)  and  rewriting 
it  as : 


Y  Y  Y 

11,  .  2  22,  .2  33,  .2 

aij  2  ('a22_a33')  2  ^a33_ail^  "  2  ^°ll"a22^ 


12 


23 


13 


=  1 


6.6 


12 


23 


23 


If  the  superscript  N  refers  to  values  found  by  nodal  equilibrium  equa¬ 
tions  and  E  refers  to  the  value  found  by  the  elements  equilibrium 
equations,  the  final  stresses  are  solved  for  by  the  following  equations: 
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o 


11 


a  =  a*  -  a*  +  aE 

22  22  33  33 


a33  a33 


6.7 


ai2  =  ° 


_N 

12 


23 


°23 


13 


aN 

13 


This  analysis  was  then  programmed  and  the  results  compared  favor¬ 
ably  with  an  experiment.  These  findings  and  conclusions  are  discussed 
in  the  next  section. 
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VII.  RESULTS  AND  CONCLUSIONS 


In  order  to  check  the  accuracy  of  the  previous  analysis,  runs  were 
made  with  this  program  that  inputted  the  requirements  of  the  following 
experiment.  The  example  used  to  check  the  program's  accuracy  was  a 
three-layered  laminated  plate. **  The  top  and  bottom  layers  were  1020 
steel,  12.7  mm  and  6.35  mm  thick,  respectively.  The  middle  layer  was 
2024  aluminum,  12.7  mm  thick.  The  plate  was  22.86  cm  by  45.72  cm.  The 
ends  were  simply  supported  and  a  large  impulsive  load  was  detonated  on 
the  top  of  the  plate.  The  results  of  the  experiment  compared  favorably 
to  another  computer  code  (acronym  HEMP) .  The.  HEMP  code  is  a  Lagrangian 
finite  difference  technique  that  utilised  an  elastic-perfectly  plastic 
model  for  solids.  It  is  relatively  time  consuming  and  costly  to  run, 
thus  the  need  for  a  more  efficient  technique  such  as  the  one  previously 
described. 


The  input  to  this  program  which  analyzed  three-dimensional  impul¬ 
sively  loaded  plates  (acronym  TIP)  is  shown  in  Appendix  A  and  is  fol¬ 
lowed  by  the  actual  input  used.  As  was  mentioned  in  Section  2,  in 
order  to  find  the  external  force,  a  special  subroutine  must  be  input 
which  gives  the  external  pressure  at  specific  locations.  For  this  ex¬ 
ample  a  list  of  data  was  given  for  the  external  pressure  at  radial 
distances  from  the  center  of  the  plate  at  discrete  time  intervals. 
Figure  8  shows  the  data  graphically  for  several  times.  Rather  than  in¬ 
put  this  data  and  extrapolate  answers,  the  pressure  was  averaged  over 
the  plate  at  discrete  time  intervals  using  the  formula: 


ave 


N 

=  l 

i=l 


,  2  2  . 
p. (r . -r.  ) 

ri  1  1-1 


LN 


7.1 


where  p^  is  the  pressure  at  radial  distance  r^. 

This  gives  an  area  type  average  and  Figure  9  shows  graphically  how 
the  results  look.  Pave  was  then  input  in  SUBROUTINE  DISFOR;  and  knowing 
the  time,  the  proper  pressure  is  extrapolated  from  these  average  values. 
This  cuts  down  considerably  on  the  necessary  input  into  SUBROUTINE 
DISFOR. 

With  the  DISFOR  and  the  input  shown  at  the  end  of  Appendix  B,  the 
program  was  executed  with  an  elastic-plastic  yield  subroutine  and  exe¬ 
cuted  again  with  an  elastic-viscoplastic  subroutine. 
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In  the  program  with  the  elastic-plastic  subroutine,  the  results 
shown  in  Figure  10  were  achieved  when  using  a  time  step  of  two  and  one 
half  microseconds  for  the  first  twenty  microseconds,  then  five  micro¬ 
seconds  for  twenty  to  sixty  microseconds,  and  ten  microseconds  for  the 
remaining  time.  There  was  a  four  by  six  grid  on  the  plate  which  is 
quite  simple.  The  results  shown  here  have  a  maximum  error  of  five  per¬ 
cent,  which  in  all  cases  was  less  than  the  error  of  the  HEMP  code. 

This  lends  credibility  to  this  method  of  analysis. 

Figure  11  gives  the  results  of  an  elastic-viscoplastic  analysis 
for  this  example.  Since  there  is  damping  in  this  type  of  analysis  the 
elastic-viscoplastic  analysis  was  stiffer  than  the  elastic-plastic 
analysis . 

Various  techniques  could  be  used  to  improve  the  error  found  in  this 
example.  Smaller  time  steps  could  be  used;  and  this  would  reduce  the 
error  but  increase  the  computation  time.  A  smaller  grid  with  more  nodal 
points  could  be  used  but  it  must  be  remembered  that  the  more  degrees  of 
freedom  a  problem  has,  the  more  mode  shapes  it  will  have  and  conse¬ 
quently  the  more  chances  of  it  going  into  an  unstable  mode  shape.  Also 
the  approximations  used  in  the  large  deflection  analysis  assume  an  ele¬ 
ment  that  has  plate-like  characteristics, 

This  is  a  highly  efficient  way  to  handle  a  very  difficult  problem. 
Its  advantages  are  in  the  speed  and  accuracy  it  produces.  The  element 
equilibrium  equations,  although  simple  in  their  calculation  of  the  ele¬ 
ment  accelerations,  allow  for  minor  error  corrections. 
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*Note:  Some  minor  subroutines  which  are  called  repeatedly  are 

omitted  for  clarity. 


Figure  1.  Flow  Diagram  of  TIP/I  Program 
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Figure  2.  Coordinates  and  Global  Numbering  System. 
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Figure  3.  Nodal  Numbering  System  for  Quadrilateral  and 
Triangular  Elements. 
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Figure  4.  Relation  Between  Global  and  Orthotropic  Properties. 
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Figure  5.  Model  Representing  the  Predicted  Stiffness  in  the 
Thickness  Direction. 
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Figure  6.  Section  of  a  Laminated  Plate  Used 
in  Additional  Element  Equilibrium 
Equations 
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Figure  7.  Definitions  of  Interlaminar  Normal  Stresses. 
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Pressure  (Mburs) 


Figure  8.  Actual  Pressure  on  Plate  at  Discrete  Time 

Intervals 
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Figure  9.  Idealized  Pressure  on  Plate  at  Discrete  Time 
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Figure  10.  Comparison  of  Results  with  Three  Time  Increments  for  the 
Elastic-Plastic  Model. 
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Deflection  in  Centimeters 


External  Pressure 


Figure  11. 


Time  in  Microseconds 


Comparison  of  Results  with  Three  Time  Increments  for  the 
Elastic-Viscoplastic  Model. 


h=  2.5  ysec 
h-  5.0  ysec 
h=  10.0  ysec 


t=  1,  20  ysec 
t=  20,  60  ysec 
t=  60,100  ysec 
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APPENDIX  A 


COMPUTER  PROGRAM  INPUT  CARD  DESCRIPTION 
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TITLE  C m. 


Format  (20a4)  Title  (Title  for  particular  case) 


CONTROL  CARD 


Format  (4 15) 

Columns  1-5 

6-10 

11-15 

16-20 


PRINT  CARD 


NUMMAT  (Number  of  different  materials;  6 
maximum) 

NUMLA  (NumDer  of  layers;  12  maximum) 

NLINC  (Number  of  load  increments  with  time; 
NLINC21 ) 

IPLOT  (Plot  parameter,  1  if  plot  required,  0  if 
no  plot  required) 


This  card  controls  the  output  that  is  generated  so  that  the 
deflections  are  printed  every  NPRINT  time  increments. 

Format  (15) 

Columns  1-5  NPRINT 


TIME  INCREMENT  CARD 


Format  (F10.5,5E10.4) 


Columns  1-10 

11-20 

21-20 

31-40 

41-50 

51-60 


BET  (3  ,  acceleration  parameter  or  Newmarx's 
parameter,  0.25) 

First  time  step  to  be  used 
Second  time  step  to  be  used 
Third  time  step  to  be  used 

Time  at  which  program  begins  using  second  time 
step 

Time  at  which  program  begins  using  third  time 
step 


MESH  GENERATION  CONTROL  CARD 


Format  (515) 

Columns  1-5 
5-10 
10-15 
16-20 
21-25 


MAXI  (Maximum  value  of  I  in  mesh;  25  maximum) 
MAXJ  (Maximum  value  of  J  in  mesh;  100  maximum) 
NSEG  (Number  of  line  segment  cards) 

NBC  (Number  of  boundary  condition  cards) 

NMTL  (Number  of  material  block  cards) 
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L1N£  SEGMENT  CARDS 


The  order  of  line  segment  cards  is  immaterial,  except  when  plots  are 
requested;  in  this  case,  the  line  segment  cards  must  define  the 
perimeter  of  solid  continuously.  The  order  of  line  segment  cards 
defining  internal  straight  lines  is  always  irrelevant. 

Format  (3(213, 2F8. 3) ,15) 


Columns  1-3 

I 

coordinate 

of 

1st 

point 

4-6 

J 

coordinate 

of 

1st 

point 

7-14 

a 

coordinate 

of 

1st 

point 

15-22 

Z 

coordinate 

of 

1st 

point 

23-25 

I 

coordinate 

of 

2nd 

point 

26-28 

J 

coordinate 

of 

2nd 

point 

29-36 

a 

coordinate 

of 

2nd 

point 

37-44 

Z 

coordinate 

of 

2nd 

point 

45-47 

I 

coordinate 

of 

3rd 

point 

48-50 

J 

coordinate 

of 

3rd 

point 

51-58 

a 

coordinate 

of 

3rd 

point 

59-66 

z 

coordinate 

of 

3rd 

point 

67-71 

Line  segment 

type  parameter 

If  the  number  in  column  71  is: 

0  Point  (input  only  1st  point) 

1  Straight  line  (input  only  1st  and  2nd  points) 

2  Straight  line  as  an  internal  diagonal  (input 
only  1st  and  2nd  points) 

3  Circular  arc  specified  by  1st  and  2nd  points  at 
the  mid-point  of  the  arc 

4  Circular  arc  specified  by  1st  and  2nd  points  at 
the  ends  of  the  arc  with  the 

coordinates  of  the  center  of  tne  arc  given  as 
tne  3rd  point  (delete  I  and  J  for  3rd  point) 

5  Straignt  line  as  a  boundary  diagonal  for  whien  1 
of  1st  point  is  minimum  for  its  row  and/or  i  of 
2nd  point  is  minimum  for  its  row  (input  only  1st 
and  2nd  points) 

6  Straight  line  as  a  boundary  diagonal  for  wnicn  1 
of  1st  point  and/or  2nd  point  is  maximum  for  its 
row  (  input  only  1st  and  2nd  points) 

Note:  In  specifying  a  circular  arc,  the  points  are 

ordered  such  that  a  counter-clocxwise  direction 
about  the  center  is  obtained  upon  moving  along 
the  boundary . 
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BOUNDARY  CONDITION  CARDS 


Each  card  assigns  a  boundary  condition  code  to  a  block  of  succes¬ 
sive  nodal  points  starting  with  N1  and  ending  with  N2,  inclusive.  (If 
code  is  0  no  card  necessary.)  Displacements  and  velocities  can  be 
prescribed  at  any  node  for  a  particular  time  through  Subroutines  BCDEL 
and  BCVEL. 

Format  (215,110) 

Columns  1-5  Starting  node  number  N1 
6-10  Ending  node  number  N2 

11-20  Boundary  condition  code 

If  the  number  in  columns  11-20  is: 

0  displacement  is  not  prescribed  (program  assigns  this  code 
automatically) 

1  displacement  is  prescribed  in  direction 

2  displacement  is  prescribed  in  X2  direction 

3  displacement  is  prescribed  in  Xg  direction 

4  displacement  is  prescribed  in  x^  and  X2  directions 

5  displacement  is  prescribed  in  X2  and  Xg  directions 

6  displacement  is  prescribed  in  Xg  and  Xj  directions 

7  displacement  is  prescribed  in  x^,  X2  and  Xg  directions 

8  through  14,  these  codes  parallel  codes  1  through  7  except 

velocities  are  prescribed  instead  of  displacements  for 
chosen  nodes 

MATERIAL  BLOCK  ASSIGNMENT  CARD 

Each  card  assigns  a  material  definition  number  to  a  block  of 
elements  defined  by  the  I,  J  coordinates.  Two  cards  for  each  layer. 


Card  1 


Format  (I5,3F10.0) 


Columns 


Card  2** 


1-5  Material  definition  number  (1  through  6) 

6-15  Material  principal  property  inclination  angle 
BETA  in  Xj  -  X£  plane 

16-25  Material  principal  property  inclination  angle 
ALPHA  in  x'^  -  x^  plane 
26-35  Bingham  viscosity  in  this  layer 


Format  (6E12.6) 


Columns 

1-12 

Yield 

stress 

in 

x' 

'l 

direction 

13-24 

Yield 

stress 

in 

X' 

1 

9 

direction 

25-36 

Yield 

stress 

in 

X' 

•I 

direction 

37-48 

Yield 

stress 

in 

X' 

-1 

x'2  direction 

49-60 

Yield 

stress 

in 

X' 

'2 

x'g  direction 

■ 

61-72 

Yield 

stress 

in 

X' 

r 

1 

x'g  direction 
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Note:  If  material  is  isotropic,  place  yield  stress  in 

first  12  spaces  and  leave  the  rest  blantc.  When 
the  program  finds  blanxs  in  the  Id-24  spaces  it 
assumes  the  material  is  isotropic  and  calculates 
the  other  yield  stresses. 

MATERIAL  PROPERTY  INFORMATION  CARDS 

Tne  following  group  of  cards  must  be  specified  for  each  material 
(Maximum  of  6). 


a.*.  MATERIAL  IDENTIFICATION  CARD 
Format  (I5,F10.0) 

Columns  1-5  Material  identification  number 

6-15  Mass  density  of  material  (if  required) 

b^_  MATERIAL  PROPERTY  CARDS 


First  Card 

Format  (6F10.0) 

Columns  1-10 
11-20 
21-30 
31-40 
41-50 
51-60 

Second  Card 


Format  (3F10.0) 

Columns  1-10 
11-20 
21-30 

LAYER  TriiCKNESS  CARD 
Format  (12F5-3) 

Columns  1-5 
6-10 
11-15 
etc . 


Modulus  of  elasticity,  E^ 
Modulus  of  elasticity,  E^ 
Modulus  of  elasticity,  E 
Poisson's  ratio, 

Poisson's  ratio,  V 
Poisson's  ratio,  V 


Shear  Modulus, 
Shear  Modulus, 
Shear  Modulus,  G^-j 


TR(1)  (ThicKness  of  layer  1) 

TH(2)  (Thickness  of  layer  2) 

TH(3)  (Thickness  of  layer  3) 

up  to  TH  (NUMLA) 
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PLOT  TITLE  CARD* 


Format  (8A10) 

Columns  1-80  Title  (Title  printed  under  each  plot) 
PLOT  GENERATION  INFORMATION  CARD* 

Format  (2F10.0) 

Columns  1-10  RMAX  (Maximum  x-^  coordinate  of  mesh) 
11-20  ZXAX  (Maximum  X2  coordinate  of  mesh) 

*Note:  use  only  if  IPLOT  =  1  (plot  required) 


INITIAL  CONDITION  CARDS 


Card  1 

Format  (15) 

Columns  1-5  INIDV  number  of  initial  displacement  and 

velocities  cards  (see  Card  2)  specifying 
initial  displacements  and  velocities 


Card  2 

Each  card  assigns  an  initial  displacement  and/or  initial  velocity 
to  a  specific  nodal  point.**  The  number  of  these  cards  is  equal  to 
INIDIV. 

Format  (I5,2E12.6) 

Columns  1-5  Nodal  Point 

6-17  Initial  Displacement 

18-29  Initial  Velocity 

**Note:  Used  only  if  there  are  initial  displacements  or  velocities, 
otherwise  the  program  initializes  these  values  as  zero  and  no  input 
is  necessary. 
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APPENDIX  B 


PROGRAM  TIP  FOLLOWED  BY  ITS  INPUT 
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PROGRAM  TIPI  (INPUT,  OUTPUT,  DATP ,TAP£5=INPUT , TAPE6=OUTPUT ,TAPE1  , 
1TAPE2,TAPE3,TAPE4=DATP) 

C*  »*»»»»»»««****«««*««**««*««*»****«* 

C  THREE  DIMENSIONAL*’ ANALYSIS  OF  AN  IMPULSIVELY  LOADED  LAMINATED  PLATE 
C  BY  THE  FINITE  ELEMENT  METHOD  WITH  ORTHOTROPIC  PLASTIC  YIELDING 
C  AND  NONLINEAR  STRAINS 

Q#  *********************************** 

INTEGER  CODE 

COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /NPR/NPRINT 

COMMON /MATP/RO(12) ,E(9, 12) ,EE(9) ,ETA( 12) 

COMMON /ARG/XXXC 10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 , 6 ) , XI (10) ,SIG( 12)  ,N  ,M 
COMMON /NPDATA/XC 200) ,Y(200) ,C0D£(200) ,NPNUM( 10 ,20) 

COMMON /ELDATA/BETA( 12) ,ALPHA(12) ,TH(12) ,IX(200 ,4) ,MATRIL( 12) 

COMMON /RESULT/D (6, 6) ,C(6 , 6 ) ,CNS (6 , 6 ) 

COMMON /TD/IMIN (100) ,IMAX( 100 ) , JMIN (25 ) , JMAX(25 ) ,MAXI ,MAXJ ,NMTL ,NBC 
COMMON /BASI C2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM, TSUM( 6 ) , SIGMA ( 12,50,6) ,DSIG(6 ) , F (600 ) 

COMMON /DISP/DELN 1 (600 ) , DELN (600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/A( 600) ,B(600 ) ,CM(8 ) 

COMMON /FIRS T/FO( 600 ) ,DER( 600) 

COMMON /DELT /XDEL ( 4 , 1 8 ) 

COMMON /DELTRI /DELTA ( 4 , 1 8 ) 

COMMON /PLY LD / SIGY (12,6) ,DEPS(6 ) ,EPS( 12 , 50 , 6 ) 

DIMENSION  TITLE (20) 

0*  it****#****#**************#**####### 

C  READ  AND  WRITE  CONTROL  INFORMATION 

C*  *********************************** 

50  READ(5 , 1000) TITLE, NUMMAT, NUMLA, NLINC,IPLUT 
IF(EOF(5 ) )920 ,88 
88  READ(5,1004)NPRINI 

IF ( EOF ( 5 ) . NE . 0 . ) GO  TO  920 

WRITE ( 6 , 2000 ) TITLE , NUMLA , NUMMAT , NLINC 

WRITE(4)NLINC, NUMLA 

READ (5, 1002)  BET, (HH(I) ,1=1 ,3) ,(HT(J) ,J=1 ,2) 

WRITE (6 , 2001 )BET ,(HH(I),I=1,3)»(HT(J),J=1 ,2) 

H=HH(1 ) 

HN=0 .0 

TIME=0.00 

NTIME=0 

C*  *********************************** 

C  GENERATE  FINITE  ELEMENT  MESH 

100  CALL  MESH 
MPRINT=0 

DO  230  N=1, NUMNP 
IF(MPRINT.NE.O)  GO  TO  220 
WRITE(6 , 2003 ) 

MPRINT=59 
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220  MPRINT=MPRINT-1 

230  WRITE(6 ,2004)  N,X(N),Y(N) 

IT=NUMNP*(NUMLA+1)*3 
440  MPRINT=0 

DO  460  N=1 , NUMEL 
IF(MPRINT.NE.O)  GO  TO  450 
WRITE (6, 2008) 

MPRINT=59 

450  MPRINT=MPRINT-1 
IX =IX(N , 1 ) 

JJ=IX(N,2) 

KK=IX(N , 3) 

LL=IX(N ,4) 

460  WRITE(6 ,2009 )  N, (IX(N,I) ,1=1 ,4) 

C*  *********************************** 

C  READ  AND  WRITE  MATERIAL  PROPERTIES 
c*  *********************************** 

500  CONTINUE 

DO  510  M=1 ,  NUMMAT 

READ(5 , 1004)  MTYPE , (RO(MTYPE) ) 

WRITE (6 ,2010)  MTYPE,  RO( MTYPE) 

READ(5, 1005 )(E(J, MTYPE ),J=1 ,9) 

WRITE(6,2011)(E(J, MTYPE) ,J  =  1  ,9) 

510  CONTINUE 

READ (5 , 1006 ) (TH(I ) ,1=1 , NUMLA) 

WRITE (6 , 1007  ) 

1007  FORMAT ("OTHICKNESSES") 

WRITE (6 , 1006 ) (TH(I ) ,1=1 , NUMLA) 

WRITE(6 , 1008 ) 

1008  FORMAT ("OYIELD  STRESSES") 

WRITE (6 , 1009 ) ( (SIGY(I , J) ,J=1 ,6) ,1=1 , NUMLA) 

1009  FORMAT(6 (2X ,E15 .7) ) 

WRITE (6 ,1108) (ETA(I) ,1=1 , NUMLA) 

1108  FORMAT (//"  COEFFICIENTS  OF  VISCOSITY"/6(2X,E12.6)/) 

CALL  USTART 

CALL  UDIMEN(10. , 10. ) 

IF  (IPLOT.EQ. 1 )  CALL  MPLOT 
DO  800  1=1, NUMLA 
DO  800  J=1 ,NUMEL 
DO  800  K=1 ,6 
SIGMA (I,J,K)=0.00 
800  CONTINUE 
CALL  INIT 
DO  900  NL=1 ,NLINC 
IF(NL.GT.I)  GO  TO  721 
DO  720  N=1, NUMLA 
ALPHA(N)=ALPHA(N)/57. 295780 

720  BETA(N)=BETA(N)/57. 295780 

721  CONTINUE 
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C  FORM  STIFFNESS  MATRIX 
C 

DO  850  1=1,4 
DO  850  J=1 , 18  ' 

DELTA(I, J)=0.00 
850  CONTINUE 
CALL  DIFF 
900  CONTINUE 
CALL  UEND 
910  GO  TO  50 

1000  FORMAT (20A4/6I5,F5. 0,515) 

1001  FORMAT(3F10.0) 

1002  FORMAT (F10.5,5E10.4) 

1004  FORMAT  (I5,F10.0) 

1005  F0RMATC6F10.0) 

1006  FORMAT  (12F5.3) 

2000  FORMAT  (2H1  ,20A4/ 

1  33H0  NUMBER  OF  LAYERS - 14/ 

2  33H0  NUMBER  OF  MATERIALS - 14/ 

3  33H0  NUMBER  OF  LOAD  INCREMENTS - 14/) 


2001  F0RMAT(4lH0  ACCELERATION  PARAMETER,  BETA - FI 0.5/ 

148H0  TIME-STEP  SIZE.HH - 3 (2X, El 0.4)/ 

251  HO  TIME  TO  CHANGE  TIME  STEP  SIZE - 2 (2X, El 0.4) 

3  /) 

2003  FORMAT  (35H1  N  X  Y  ) 

2004  FORMAT  (I5.2F10.4) 


2008  FORMAT  (51H1  EL  I  J  K  L  ANGLE  BETA  ANGLE  ALPHA) 

2009  FORMAT  (I5,4l4,2F13.3) 

2010  FORMAT  ( 1 HI , "MATERIAL  IDENTIFICATION  NUMBER  =",12/ 

21H  , "MASS  DENSITY  =",E15.7) 

2011  FORMAT  ( 

11H  , "MODULUS  OF  ELASTICITY-EN  =",E15.7/ 

21H  , "MODULUS  OF  ELASTICITY-ES  =",E15.7/ 

31H  ."MODULUS  OF  ELASTICITY-ET  =",E15.7/ 

41H  ."POISSON  RATIO-NUNS  =",E15.7/ 

51H  ."POISSON  RATIO-NUNT  =",E15.7/ 

61H  ."POISSON  RATIO-NUST  =",E15.7/ 

71H  ."SHEAR  MODULUS-GNS  =",E15.7/ 

8lH  ."SHEAR  MODULUS-GST  =",E15.7/ 

91H  ."SHEAR  MODULUS-GTN  =",E15.7/) 

2016  FORMAT  (26H  THE  SYSTEM  CONVERGED  IN  I2.11H  ITERATIONS) 

2017  FORMAT  (33H  THE  SYSTEM  DID  NOT  CONVERGE  IN  I2.11H  ITERATIONS) 
920  STOP 

END 

SUBROUTINE  ANGLE  (R.Z.RC.ZC.ANG) 

INTEGER  CODE 

COMMON/BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/MATP/RO(12),E(9,12),EE(9),ETA(12) 

COMMON /ARG /XXX ( 10 ) ,YYY (10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 ,6) ,XI(10) ,SIG(12) ,N,M 
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C0MM0N/NPDATA/X(200) ,Y(200 ), CODE (200 ) ,NPNUM( 10 , 20 ) 

COMMON /ELD ATA/BETA( 12 ) , ALPHA ( 12) ,TH(12) ,  IX (200 , 4 ) ,MATRIL( 12 ) 
COMMON/RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

COMMON /  TD  / 1  MIN,(  1 0  0 ) ,IMAX(100) ,JMIN(25) , JMAX(25 ) ,MAXI ,MAXJ ,NMTL ,NBC 
COMMON /BASIC2/ BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM , TSUM ( 6 ) , SIGMA (12,50,5), DSIG ( 6 ) , F ( 6  00 ) 
COMMON/DISP/DELN 1 (600) ,DELN(600) , DEL (600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/A( 600) ,B(600) ,CM(8) 

c#  *********************************** 

C  FIND  ANGLE  OF  INCLINATION  BETWEEN  0  AND  2*PI 

c#  *********************************** 

PI=3. 1415927 
D1 = (Z-ZC) 

D2= (R-RC) 

IF(ABS(R-RC) .GT.1 .E-8)  GO  TO  100 
ANG=PI/2. 

IF(D1.GT.1.E-8)  RETURN 

ANG=-ANG 

RETURN 

c*  *********************************** 

C  ALLOW  CIRCLE  TO  CROSS  AXIS 

c*  ******************  ***************** 

100  ANG=ATAN2(D1 ,D2) 

RETURN 

END 

SUBROUTINE  BAND  (A,B,NROW,NCOL,NMID,NRES,KO) 

DIMENSION  A(13,3),B(13,1) 

NL=NMID-1 
NR=NCOL-NMID 
NRM=NR0W-1 
DO  2  1=1, NRM 

IF ( A ( I , NMID ) . EQ . 0 . 0 )  GO  TO  7 
RE=1 . 0/A(I ,NMID) 

LL=NROW-I 

LR=LL 

IF(LL.GT.NL)  LL=NL 
IF(LR.GT.NR)  LR=NR 
DO  2  J=1 ,LL 
JR=I+J 
JC=NMID-J 
RA=-A( JR , JC)*RE 
DO  1  K=1 ,LR 
KC=JC+K 
IC=NMID+K 

1  A( JR ,KC) =A( JR ,KC)+RA*A(I , IC) 

DO  2  K=1 ,NRHS 

2  B(JR,K)=B(JR,K)+RA*B(I,K) 

if(a(nrow,nmid) .eq.o.o)  go  to  6 

RE=1 . 0/A(NR0W,NMID) 

DO  3  1=1 ,NRHS 
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o  o  o  o 


3  B ( NRO W , I ) =B ( NROW , I )  *RE 
DO  5  1=1, NRM 
IR=NROW-I 
RE=1./A(IR,NMID) 

LR=I 

IF(LR.GT.NR)  LR=NR 
DO  4  J = 1 , LR 
JR=IR+J 
JC=NMID+J 
DO  4  K=1 ,NRHS 

4  B(IR,K)=B(IR,K)-A(IR,JC)*B(JR,K) 

DO  5  J=1 ,  NRHS 

5  B(IR,J)=B(IR,J)*RE 
K0=0 

RETURN 

6  I=NROW 

7  KO=I 

WRITE(6 , 100 )  KO 
RETURN 

100  FORMAT (//"  YOU  GOOFED  -  THERE  IS  A  ZERO  ON  THE  PRINCIPAL  DIAGO 

INAL  IN  THE", 14,"  TH  ROW."//) 

END 

SUBROUTINE  BCDEL(I) 

COMMON /BASI C2 / BET , H , HN , HH ( 3 ) , HT (2 ) , TIME , NLA Y , IF LAG , IT , NTIME 
COMMON /DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 

c  i  *  «  «  i  *  *  *  t  *  t  *  *  i  *  «  «  *  «  «  *  «  i  «  i  i  i  «  *  «  «  «  *  *  « 

IF  DISPLACEMENTS  ARE  GIVEN  AT  SPECIFIED  NODES,  THIS  IS  THE  POINT  IN 
THE  SUBROUTINE  WHERE  THAT  SPECIFIC  INPUT  IS  GIVEN  WITH  RESPECT  TO 
TIME  AND  NODAL  LOCATION 

*»tttttt*tt***«********  ****######**### 

DEL (I) =0.0 
RETURN 
END 

SUBROUTINE  BCVEL(I) 

COMMON / BASI C2 /BET , H , HN , HH (3 ) , HT (2 ) , TIME , NLAY , IF LAG , IT , NTIME 
COMMON/DISP/DELN 1(600) ,DELN(600) ,DEL(600) ,GNM1(600) ,GNM2(600) 
COMMON /FIRST/FO( 600 ) ,DER( 600) 

COMMON /MASS 1 /XMINV( 600 ) ,EEKM(13) ,F1 (600) ,F2(600) 

DERN=DER(I) 

C  «  •  «  «  •  *  «  •  *  «  *  «  «  «  «  i  «  i  «  «  •  «  «  «  «  i  i  *  i  •  *  (  •  *  * 

C  IF  VELOCITIES  ARE  GIVEN  AT  SPECIFIED  NODES,  THIS  IS  THE  POINT  IN  THE 
C  SUBROUTINE  WHERE  THAT  SPECIFIC  INPUT  IS  GIVEN  WITH  RESPECT  TO  TIME 
C  AND  NODAL  LOCATION 

c  i  i  *  ««  i  ««  i  ««*««*«»*«*«*««*  i  i  *  *  «  *  *  *  «  * 

DER(I)=0.0 

DEL(I) =DELN ( I )+( 1 . -2 . *BET ) *H*DERN+2 . *BET*H*DER(I)+( . 5-2 . *BET ) 

1  *H**2*XMINV(I)«(GNM1(I)+F1(I)) 

RETURN 

END 

SUBROUTINE  CIRCLE ( ANG 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 
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INTEGER  CODE 

COMMON /BASIC/VOL ,  NUMNP , NUMEL , NUMLA , NCG 
COMMON /MATP/RO( 12) ,E(9,12) , EE(9 ) ,ETA( 12 ) 

COMMON /ARG/XXX,(  10)  ,YYY(10)  ,S(24)  ,XX(3)  ,YY(3) , 

1CRZ(6 ,6) , XI (10)  ,SIG(12) ,N,M 

C0MM0N/NPDATA/X(200 ) ,Y(200) ,C0DE(200) ,NPNUM( 10 ,20) 

COMMON /ELD AT A / BET A ( 1 2 ) ,ALPHA(12) ,TH( 12) ,IX(200 ,4) ,MATRIL(12) 

COMMON /RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

COMMON /TD/ I MIN (100) ,IMAX(100) , JMIN(25) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT (2 ) , TIME , NLAY , IFLAG , IT , NTIME 
C0MM0N/SIGM/BSUM,TSUM(6) ,SIGMA( 12,50,6) ,DSIG(6) ,F(600) 
C0MM0N/DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
C0MM0N/MASS/A(600) ,B(600) ,CM(8) 

DIMENSION  AR(10,  40),AZ(10,  40) 

EQUIVALENCE  (X ( 1 ) , AR ),(Y(1), AZ ) 

q*  *********************************** 

C  FIND  INTERSECTION  OF  LINE  AND  CIRCLE  =  NEW  R  AND  Z 

C*  *********************************** 

ANG1=ANG1+DELPHI 

RR=SQRT ( (RSTRT-RC) **2+( ZSTRT-ZC) **2 ) 

AR ( I , J ) =RC+RR*COS ( ANG 1 ) 

AZ(I , J)=ZC+RR*SIN(ANG1 ) 

RETURN 

END 

SUBROUTINE  DIFF 
INTEGER  CODE 

COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/MATP/RO( 12 ) ,E(9 , 12) ,E£(9 ) ,ETA( 12 ) 

COMMON /ARG/XXX( 10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 ,6) ,XI(10) ,SIG(12) ,N,M 

COMMON/NPDATA/X(2O0) ,Y(200) ,C0DE(200) ,NPNUM( 10 ,20) 
COMMON/ELDATA/BETA( 12) ,ALPHA(12) ,TH(12) , IX (200 , 4 ) ,MATRIL( 12 ) 

COMMON /RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

C0MM0N/TD/IMIN(100) ,IMAX(100) ,JMIN(25) ,JMAX(25) , MAXI , MAX J , NMTL , NBC 
COMMON /BASI C2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLA Y , IFLAG , IT , NTIME 
C0MM0N/SIGM/BSUM,TSUM(6) , SIGMA ( 12,50,6) ,DSIG(6) ,F(600) 

COMMON /DISP/DELN 1(600), DELN ( 600 ) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
C0MM0N/MASS/A(600) ,B(600) ,CM(8) 

C0MM0N/FIRST/F0(600) ,DER(600) 

C0MM0N/MASS1/XMINV(600) ,EEKM( 13) , FI (600) ,F2(600) 

COMMON /NPR/NPRINT 

COMMON /PL YLD/SIGY (12,6) ,DEPS(6 ) ,EPS( 12 ,50 , 6 ) 

HN=H 

TIME=TIME+H 
NTIME=NTIME+1 
DO  20  1  =  1  ,2 

IF(TIME.LT.HT(I) .0R.H.EQ.HH(I+1 ))G0  TO  20 
IF(TIME.GT.HT(2) .AND.H.EQ.HH(3))G0  TO  30 
H=HH(I+1 ) 

TIME=TIME-HN+H 
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on  o  o 


WRITE (6, 10)HN,H,TIME,HT(I) 

10  FORMATO'O  TIMES  "  ,M2X,E10.4)/) 

GO  TO  30 
20  CONTINUE 

30  IF(NTIME.GT. 1 )G0  TO  280 
WRITE (6, 999) 

999  FORMATC'I INITIAL  FORCE  VECTOR"/) 

WRITE (6 , 1000) (FO (I) ,1=1 , IT) 

WRITE (6, 998) 

998  FORMAT ("OINITIAL  DISPLACEMENTS") 

WRITE (6 , 1000) (DEL (I) ,1=1 ,IT) 

WRITE (6, 997) 

997  FORMAT ("OINITIAL  VELOCITIES") 

WRITE (6 , 1000) (DER(I) ,1=1 ,IT) 

WRITE (6, 996) 

996  FORMATC'I") 

1000  FORMAT(9(2X,E12.6)) 

280  CONTINUE 

DO  390  1=1, IT 
IF(NTIME.EQ. 1 )GO  TO  390 
F2(I)=F1 (I) 

390  FI (I )=F(I ) 

CALL  LOAD 
CALL  STRESS 
CALL  STIFF 
CALL  INT 

IF(NTIME.EQ.I)  GO  TO  98 

IF( (NTIME/NPRINT)*NPRINT.NE.NTIME)  GO  TO  99 
98  CONTINUE 

WRITE(6 ,26 ) ( (I , J , (EPS(I , J ,K) ,K=1 ,6), 1=1 ,NUMLA),J=1 ,NUMEL) 

WRITE (6 ,25 ) ( (I , J, (SIGMA(I , J ,K) ,K=1 ,6), 1=1 ,NUMLA),J=1 ,NUM£L) 

WRITE ( 6 , 1 00  3 ) NTIME , TIME 

1003  FORMAT (/"ODEFLECTION  FOR  CYCLE=",I3,"  TIME  =",E12.6/ 

1  "  NODE", 15X, "DEL  X",15X,"DEL  Y",15X,"D£L  Z"/) 

NUP=NUMLA+1 
DO  1026  K=1 ,NUP 
DO  1026  1=1 ,NUMNP 
NP=I+(K-1 )*NUMNP 
NP3=3*NP 

WRITE(6, 1001 )NP, DEL (NP3-2),DEL(NP3-2),DEL(NP3) 

1001  FORMAT (16, 3E20. 10) 

1026  CONTINUE 

MI=NUMNP*3 

WRITE (6 , 15 ) (MX, DEL (MX) ,DELN(MX) ,DELN1 (MX) ,B(MX) ,GNM1 (MX) , 

1  GNM2(MX),F(MX),MX=3,IT,MI) 

15  FORMAT ( / "  #",5X,"DEL  N",8X,"DELN  N-1 " ,5X , "DELN1  N-2",8X, 

1  "B  N" ,8X,"GNM1  N-1 " ,8X , "GNM2  N-2",8X,"F" 

2  /(I4,2X,E12.6,2X,E12.6,2X,E12.6,2X,E12.6,2X,E12.6,2X,E12.6, 

3  2X.E12.6)) 

25  FORMAT (//"I  LAYER" ,2X , "ELEMENT" ,5X , "SIGMA  X  ",6X, "SIGMA  Y  ",6X, 
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1  "SIGMA  Z  ",6X, "SIGMA  XY"  ,  6X , "SIGMA  YZ" , 6X , "SIGMA  XZ"/(2l8,2X>6El4 

2  .6)) 

26  FORMATC//"  LAYER", 2X, "ELEMENT", 5X,"  EPS  X  ",6X,"  EPS  Y  ",6X, 

1  "  EPS  Z  " , 6X ,,"  EPS  XY " , 6X , "  EPS  YZ",6X,"  EPS  XZ"/(2I8,2X,6E14 

2  .6)) 

99  CONTINUE 

C  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  'X1  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  'X1  T  'I1  T  T  T  T  T  T  T 
NTOP=NUMLA*3*NUMNP+3 
ETIM=TIME»1 .0E+06 
NL1 =3*NUMNP 
NL2= (NUMLA+1 )*NUMNP«3 

WRITE ( 4  )DEL ( 3 ) , DEL (NTOP ) , ( DEL ( I J ) , I J=4 , NL2 , NL 1 ) , ETIM 
C  ^ iji iji iji iji iji iji iji iji iji  iji  iji iji  iji <ji »ji «ji jji «ji T TT  iji iji  TT  T T T T T T T T  T T T  T T T T T T  T T T T 

IF( (NTIME/5 )*5 .NE ,NTIME)GO  TO  1039 
CALL  MPLOT 
1039  CONTINUE 
RETURN 
END 

SUBROUTINE  DISFOR 

COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /FOR /FF ( 5 ) , FC ( 5 ) 

COMMON/ARG/XXX(10),YYY(10),S(24),XX(3),YY(3), 

1CRZ(6 ,6) , XI ( 1 0 ) , SIG ( 1 2 )  , N , M 
COMMON/NUMB/NITER , HD AT , TO 
DIMENSION  DIST(5) 

DIMENSION  FF1 (19) 

DATA  (FF1 (I) ,1=1 ,19)/1 .4628564E-1 , 8.561 575E-2 ,2.231 6769E-2 , 

1  6 . 01 2404 IE-3 , 2 . 2297704E-3 ,1.1 280587E-3 , 6 . 0944 1 42E-4 , 

1  3 . 5506404E-4 , 2 . 266557  IE-4 , 1 . 5336 1 39E-4 ,1  .073691 8E-4 , 

1  8 . 1 0576 1 1E-5 , 6 . 5648474E-5 , 5 . 5056944E-5 , 5 . 0003866E-5 , 

1  4 .4950788E-5 , 3 . 9897700E-5 , 3 • 6072003E-5 , 3 • 2246295E-5/ 

11  NITER=1 9 

TNITER=42.0E-6 
HDAT=2.0E-6 
T0=6 .OE-6 


DO  10  1=1 ,5 

DIST(I)=SQRT(XXX(l)**2+YYY(I )**2) 

10  CONTINUE 

YBLAST=.2362*(TIME*1 .0E+06) 

IF (TIME. GT. TO)  GO  TO  40 
DO  45  1=1,4 

45  FF(I)=-FF1 (1 )*14.5E6 
GO  TO  18 

40  IF(TIME-TNITER)  15,16,16 

15  CALL  POINT(TIME ,NUM,TN ) 

DO  12  1=1  ,4 

FF(I)=-(FF1 (NUM)+(FF1 (NUM+1 )-FF1 (NUM) )*(TIME-TN)/HDAT)*14.5E6 
12  CONTINUE 
GO  TO  18 

16  DO  17  1=1,4 
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FF ( I ) =— (FF 1 (NITER)+(FF1 (NITER)-FFI (NITER-1 ) )*(TIME-TNITER)/HDAT) 

1  *14.5E6 

IF(FF(I).QT.0.0)  FF(I)=0.0 

17  CONTINUE 

18  CONTINUE 

FF (5 )= (FF ( 1 )+FF(2)+FF(3)+FF(4) )/4.0 
DO  30  1*1,5 

IF ( DIST ( I ) . GT . YBLAST )  FF(I)=0.0 
30  CONTINUE 
RETURN 
END 

SUBROUTINE  INIT 
INTEGER  CODE 

COMMON /MASS 1 /XMINV( 600 ) ,EEKM(13) ,F1 (600 )  ,F2(600 ) 

COMMON /FIRST/FO( 600) ,DER(600) 

COMMON/BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
C0MM0N/MATP/R0(12) ,E(9 , 12) ,EE(9) ,ETA(12) 

C0MM0N/ARG/XXX(10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6,6) ,XI(10) ,SIG(12) ,N,M 
COMMON /NPDATA/X (200), Y( 200), CODE ( 200 ) , NPNUM( 1 0 , 20 ) 

COMMON /ELDATA/BETA (12), ALPHA ( 1 2 ) , TH ( 1 2 ) , IX ( 200 , 4 ) ,MATRIL (12) 
COMMON/RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

COMMON /TD /IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT (2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON/SIGM/BSUM , TSUM ( 6 ) , SIGMA (12,50,6), DSIG ( 6 ) , F ( 600 ) 

COMMON /DISP /DELN 1(600), DELN ( 600 ) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/ A ( 600 ) ,B( 600 ) , CM( 8 ) 

COMMON /PLYLD/SIGY (12,6), DEPS ( 6 ) , EPS ( 1 2 , 50 , 6 ) 

C  LET  INITIAL  DEFLECTION,  VELOCITY,  AND  STRESS  BE  ZERO 
C  LET  INITIAL  FORCE  BE  LINC 
DO  100  1=1, IT 
DELN (I) =0.00 
DER(I)=0.00 
DEL(I)=0.0 
FI (I)=0.0 
F2(I)=0.0 
DELN1(I)=0.0 
F(I)=0.0 
B(I)=0.0 
GNM1 (I)=0.00 
GNM2(I)=0.0 
FO(I  )=0 .00 
A(I )=0 .0 
100  CONTINUE 

DO  110  1=1, NUMLA 
DO  110  J=1, NUMEL 
DO  110  K=1 ,6 
110  EPS(I , J ,K)=0 .0 
READ(5,200)INIDV 
IF(E0F(5).NE.0.0)G0  TO  220 
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200  FORMAT (215) 

READ (5 ,210 ) ( (ND ,DEL(ND) ,DER(ND) ) ,1=1 ,INIDV) 

210  FORMAT (I5,2E12.6) 

220  CONTINUE 
RETURN 
END 

SUBROUTINE  INT 
INTEGER  CODE 

C0MM0N/FIRST/F0(600 ) ,DER(600) 

COMMON/MA3S1/XMINV(600) ,EEKM(1j) ,F1 (600) ,F2(600) 

COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
C0MM0N/MATP/R0(12) ,E(9, 12) ,EE(9) ,ETA( 12) 

COMMON /ARG/XXX( 10) ,YYY(10) ,S(24) ,XX(3),YY(3) , 

1 CRZ ( 6 , 6 ) , XI ( 1 0 ) , SIG ( 1 2 ) , N , M 
COMMON /NPDATA/X(200 ) ,Y(200) ,CODE (200) ,NPNUM( 10 ,20) 

COMMON /ELD ATA/BETA( 12) ,ALPHA(12) , TH ( 1 2 ) , IX (200 ,4) ,MATRIL( 12) 

COMMON /TD/IMIN (100) , IMAX( 100 ) , JMIN (25 ) , JMAX(25 ) ,MAXI ,MAXJ ,NMTL ,NBC 
COMMON /BASIC2/ BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM, TSUM( 6 ) , SIGMA ( 12,50,6) ,DSIG (6 )  ,F(600 ) 
COMMON/DISP/DELN 1(600) ,DELN(600 ) , DEL (600 ) ,GNM1 (600 ) ,GNM2(600 ) 
C0MM0N/MASS/a(600) ,B(600) ,CM(8) 

COMMON /DISP 1 /TF ( 3  9 » 4 ) ,TFI(4,4) 

DIMENSION  THT(25) ,A0T(25) , BOT ( 4 ) ,  FOT ( 4 ) , DERT ( 4 ) , DELNT ( 4 ) ,AU(13,3)  , 
1  AFT(25,4) ,GNM1T(4) ,GNM2T(4) ,SM(4) ,DELN1T(4) ,DELT(4) ,BU(13,1) 

1  ,EKM(13) 

BH=BET  *H**2 
DIF  1=0. 50-BET 
BHN=BET*H*HN 

q  ip  /p  <p  ip  ip  >p  ip  ip  <p  ip  fji  ip  fji  ip  fp  ip  ip  ip  ip  ip  ip  <p  ip  ip  ip  »p  ip  rp  ip  fp  ip  rp  ip  i pi  ?p  ip  fp  rp  ip  ip  ip  rp  rp  fp  rp  rp  ip  ip  tp  fp  pi  pi  pi  pi  pi  pi  pi  <p  pi  pi  pi  pi  pi 

IF(NTIM£.GT.1)G0  TO  789 
EMIT=0 .0 
DHET=0 .0 

NT0P=NUMLA*3*NUMNP+3 
NL1 =NUMNP*3 
NL2=NL 1  * (NUMLA+1 ) 

WRITE(4 )DEL(3 ) ,DEL(NT0P) , (DEL(IJ ) ,IJ=4 ,NL2 ,NL 1 ) .EMIT 
789  CONTINUE 

Q  j  £  P  T  ^  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  X  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T  T 
C  H  AND  BETA  ARE  INPUT  VARIABLES  STORED  IN  COMMON 
NUP=NUMLA+1 
NPT= (NUMLA+1 )*NUMNP 
DO  100  1=1, IT 
DELN 1(1) =DELN ( I ) 

DELN(I)=DEL(I) 

100  CONTINUE 

IF ( NTIME .GT . 1 )G0  TO  180 
C  OBTAIN  INVERSE  OF  MASS  MATRIX, XMINV 
DO  120  1=1, IT 
XMINV (I)=1/A(I) 

120  CONTINUE 
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C  STARTING  PROCEDURE 

C  DELN  AND  DER  ARE  INITIAL  CONDITIONS  ON  DISPLACEMENT 
C  GNM1  IS  INPUT  VECTOR 
DO  160  K=1 ,NPT 
ID=CODE(K) 

I=3*K 

IF(ID.EQ.O)GO  TO  140 

IF ( ID . EQ . 3 • OR . ID . EQ . 5 . OR . ID . EQ . 6 . OR . ID . EQ . 7 )  GO  TO  125 
IF ( ID . EQ . 1 0 . OR . ID . EQ . 1 2 . OR . ID . EQ . 1 3 • OR . ID . EQ . 1 4 )  GO  TO  130 
GO  TO  140 
125  CALL  BCDEL(I) 

GO  TO  160 
130  CALL  BCVEL(I) 

GO  TO  160 

140  DEL ( I ) =DELN ( I ) +H*DER ( I ) 

DEL ( I ) =DEL ( I ) +XMINV ( I ) » ( BH»B ( I ) -DIF 1 *H*»2  *GNM 1(1)) 
DEL(I)=DEL(I)+XMINV(I)»(DIF1«H«2»F0(I)+BH*F(I)) 

160  CONTINUE 
GO  TO  240 
180  CONTINUE 
KK=0 

C  CALCULATE  DEL  BY  FINITE  DIFFERENCE  EQUATIONS 
DO  220  K=1 ,NPT 
ID=CODE(K) 

I=3*K 

IF(ID.EQ.O)  GO  TO  200 

IF ( ID . EQ . 3 . OR . ID . EQ . 5 . OR . ID . EQ . 6 . OR . ID . EQ . 7 )  GO  TO  185 
IF ( ID . EQ . 1 0 . OR . ID . EQ . 1 2 . OR . ID . EQ . 1 3 • OR . ID . EQ . 1 4 )  GO  TO  190 
GO  TO  200 
185  CALL  BCDEL(I) 

GO  TO  220 
190  CALL  BCVEL(I) 

GO  TO  220 

200  DEL ( I ) =DELN ( I )+ (DELN ( I )-DELN 1(1)) «H/HN 

DEL(I)=DEL(I)+XMINV(I)»(BH«B(I)+(H+HN)»H»DIF1»B(I)+BHN»GNM1(I)) 
DEL(I)=DEL(I)+XMINV(I)»(BH«F(I)+(H+HN)»H«DIF1»F1(I)+BHN»F2(I)) 
220  CONTINUE 
240  DO  480  K=1,NUMNP 
C  FIND  PREDICTED  STIFFNESS  (EEKM) 

260  IF(NTIME.GT.I)  GO  TO  300 
IF(K.GT.I)  GO  TO  300 
I 1 =IX( 1,1) 

I2=IX( 1,2) 

I3=IX( 1 ,3) 

I4=IX(1 ,4) 

AR=SQRT ( ( X ( 12 ) —X (II ) )**2+(Y(I2)-Y(I1 ) )**2) 

BR=SQRT( (X(I2)-X(I3) )**2+(Y(I2)-Y(I3) )**2) 
CR=SQRT((X(I4)-X(I3))**2+(Y(I4)-Y(I3))««2) 

FR=SQRT ( ( X ( I 4 ) —X (II ) )**2+(Y(l4)-Y (II ) )**2) 

PR=SQRT ( ( X ( I 3 ) —X (II) )**2+(Y(I3)-Y(I1 ) )**2) 
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QR=SQRT((X(I2)-X(I4))*»2+(Y(I2)-Y(I4))**2) 
AREA=SQRT(4.0»PR»*2«QR**2-(BR**2+FR**2-AR**2-CR**2)**2)/4.0 
REWIND  1 

READ(1)((CRZ(KJ,J) ,J=1 ,6) ,KJ=1 ,6) 

DO  280  1= 1 , NO ML A 

EEKM(l)=CRZ(3,3)*(AREA/4.0)/TH(l) 

DO  280  L=1,NUMEL 

280  READ ( 1 ) ( ( CRZ ( K J , J ) ,J=1 ,6),KJ=1,6) 

300  DO  340  1=1 , NUMLA 

EKM(I)=EEKM(I)*A(K«3)/A(3) 

340  CONTINUE 

C  FIND  NEW  COEFFICIENTS  DUE  TO  STIFFNESS  (KM) 

C  THIS  IS  DONE  BY  SOLVING  THE  SIMULTANEOUS  EQUATIONS  THROUGH  THE 
C  THICKNESS  BY  USE  OF  A  N0N3YMMETRIC  BANDED  MATRIX 
DO  400  1=1, NUP 
II=K*3+(I-1 )*NUMNP*3 
ID=CODE(K+(I-1 )*NUMNP) 

IF(ID.EQ.O)GO  TO  370 

IF(ID .EQ. 3 • OR. ID .EQ. 5 .OR . ID.EQ. 6 .OR . ID .EQ.7 )GO  TO  350 
IF ( ID . EQ . 1 0 . OR . ID . EQ . 1 2 . OR . ID . EQ . 1 3 • OR . ID . EQ . 1 4 )GO  TO  360 
GO  TO  370 
350  CALL  BCDEL(II) 

GO  TO  380 
360  CALL  BCVEL(II ) 

GO  TO  380 

370  IF(I.EQ.1)GO  TO  371 
IF(I.EQ.NUP)GO  TO  372 

AU(1 , 1 )=-BH»XMINV(II)*EKM(I-1 ) 

AU (1 , 2 )  =  1 . 0+BH*XMINV (II ) * ( EKM(I ) +EKM( 1-1 ) ) 

AU (I , 3 ) =-BH»XMlNV(II ) *EKM( I ) 

BU (I , 1 ) =DEL (II ) -BH»XMINV (II ) * ( EKM( I ) *DELN (II+NUMNP»3 ) - ( FKM (I ) 

1  +EKM(I-1 ))»DELN(II)+EKM(I-1 )*DELN(II-NUMNP*3 ) ) 

GO  TO  400 

371  AU(I, 1 )=0. 

AU (I , 2 ) = 1 . 0+BH*XMIN V ( II ) *EKM(I ) 

AU (I , 3 ) =-BH*XMlNV ( II ) *EKM(I ) 

BU (I , 1 ) =DEL ( II ) -BH*XMINV ( II ) *EKM(I ) » ( DELN ( II+NUMNP*3 ) -DELN (II ) ) 

GO  TO  400 

372  AU(I , 1 )=-BH*XMINV(II)*EKM(I-1 ) 

AU ( I , 2 ) = 1 . 0+BH*XMIN V (II ) *EKM(I- 1 ) 

AU(I,3)=0. 

BU(I, 1 )=DEL(II)+BH*XMINV(II)*EKM(I-1 )»(DELN(II)-DELN(II-3*NUMNP)) 
GO  TO  400 
380  AU (I , 1 ) =0 . 

AU(I ,2)=1 . 

AU(I,3)=0. 

BU(I, 1 )=DEL(II) 

400  CONTINUE 

CALL  BAND(AU,BU,NUP,3,2,1,KO) 

IF(KO.NE.O)GO  TO  1230 
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DO  420  1=1, NUP 
II=K*3+(I-1)*NUMNP»3 
420  DEL(II )=BU(I , 1 ) 

480  CONTINUE 
500  CONTINUE 
NVD=NUP*2 
DO  1220  K=1 ,NUMNP 
C  FIND  THE  CENTER  OF  GRAVITY 

IF(NTIME.GT.I)  GO  TO  780 
IF(K.GT.I)  GO  TO  780 
AT=A( (K-1 )*3+1 ) 

ZCGT=0 .0 
THT( 1 )=0 .0 
JL=NUMNP«3+(K-1  )*3+1 
DO  520  1=2, NUP 
THT(I)=THT(I-1 )+TH(I-1 ) 

ZCGT=ZCGT+THT(I )*A( JL) 

AT=AT+A(JL) 

520  JL=JL+NUMNP*3 
ZCG=ZCGT/ AT 
DO  560  1=1 ,NUMLA 
560  IF(ZCG.LT.THT(I+1 ))G0  TO  580 
580  NCG=I 

WRITE(6 , 600 )NCG 

600  FORMATC'  CENTER  OF  GRAVITY  IS  IN  THE  " , 15 , "  TH  LAYER") 
C  FIND  THE  TRANSFORM  MATRIX 
DO  620  1=1, NVD 
DO  620  J  =  1  ,4 
TF(I,J)=0.0 
IF  (I.GT.4)  GO  TO  620 
TFI(I , J)=0 .0 
620  CONTINUE 

DO  660  1=1, NUP 
II=2*(I-1 ) 

DO  660  J=1 ,2 
JJ=2*( J-1 )+1 
TF(II+J, JJ)=1 .0 
KK= (J-1 )*2+2 

660  TF(II+J ,KK)= (THT (I)-ZCG) 

11=1 

IF  (TF(1,2).EQ.TF(3,2))  GO  TO  700 
DO  680  1=1,2 

TFI (II  ,I)  =  1.0-(TF(1  ,2)/(TF(1 ,2)-TF(3 ,2) ) ) 

J=II+1 

TFI(J,I)=-1 .0/(TF(3,2) -TF  ( 1,2)) 

111=1+2 

TFI (II  ,III)=-TF( 1 ,2)/(TF(3 ,2)-TF( 1,2)) 

11=11+2 

680  TFI( J , III )  =  1 .0/(TF(3,2) -TF (1,2)) 

GO  TO  760 
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700  DO  720  1=1,2 

TFI(II,I)=1.0-(TF(5,2)/(TF(5,2)-TF(7,2))) 

J=II+1 

TFI(J,I)=-1 . 0 / ( TF ( 7 , 2 ) -TF (5,2)) 

111=1+2 

TFI(II ,III)=-TF(5 , 2)/(TF(7 ,2)-TF(5 ,2) ) 

11=11+2 

720  TFI(J,III)=1 .0/(TF(5 ,2)-TF(7 ,2) ) 

760  CONTINUE 

C  TRANSFORM  FORCES  AND  INITIAL  DISPLACEMENTS  AND  VELOCITIES 
780  DO  800  1=1,4 
BOT(I)=0.0 
FOT(I)=0.0 
DELN1T(I)=0.0 
DERT(I)=0.0 
DELNT(I)=0.0 
GNM1T(I)=0 .0 
800  GNM2T(I)=0.0 
DO  900  J=1  ,4 
1=0 

DO  880  IJ  =  1  ,NUP 
DO  880  L=1 ,2 
1=1+1 

N J=L+j*(K-1 )+(IJ-1 )*NUMNP*3 
IF  (J.GT.1)  GO  TO  820 
AOT(l)=A(NJ) 

820  BOT ( J )=B(NJ )*TF(I , J)+BOT( J ) 

FOT ( J )=FO(NJ )*TF(I , J)+FOT ( J ) 

GNM1T(J)=GNM1 (NJ)*TF(I, J)+GNM1T(J) 

IF  (NTIME.LT. 2)  GO  TO  860 

GNM2T ( J ) =GNM2 ( N J ) *TF ( I , J ) +GNM2T ( J ) 

860  CONTINUE 
880  CONTINUE 
900  CONTINUE 

DO  960  J=1 ,4 
1=0 

DO  940  IJ=1 ,2 
DO  940  L=1 ,2 
1=1+1 

NJ  =  (K-1  )*3+L+(IJ-1  )*NUMNP*3 
DERT( J )=DEfl(NJ )*TFI( J , I)+D£RT ( J ) 

DELNT ( J ) =DELN (N J )*TFI ( J , I)+DELNT ( J ) 

IF  (NTIME.LT. 2)  GO  TO  920 

DELN1T ( J ) =DELN 1 (NJ )*TFI( J , I)+DELN 1T( J ) 

920  CONTINUE 
940  CONTINUE 
960  CONTINUE 
C  FIND  NEW  MASS  MATRIX 
DO  980  1=1, NVD 
DO  980  J  =  1  ,4 
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980  AFT ( I , J ) = AOT ( I ) *TF ( I  ,J) 

DO  1000  1=1,4 
1000  SM(I)=0.0 

DO  1120  1=1,4 
DO  1120  J=1,NVt> 

1120  SM(I)=SM(I)+TF( J , I)*AFT(J ,1) 

C  FIND  DELT 

BH=BET*H**2 

DIF=H**2*( 1 ,0-2.0*BET) 

DIF 1=. 5-BET 
DO  1180  1=1,4 
IF(NTIME.GT.I)  GO  TO  1140 

DELT ( I ) =DELNT ( I ) +H*DERT ( I ) + ( BH*BOT ( I ) -DIF 1 *H**2*GNM1T ( I ) +DIF 1 
•  1  *H**2»F0T(I))/3M(I) 

GO  TO  1180 

1 140  DELT(I )=DELNT(I )+(DELNT (I)-DELNIT(I ) )*H/HN 

DELT ( I ) =DELT ( I ) + ( BH*BOT ( I ) + ( H+HN ) *H*DIF 1 *GNM 1 T ( I ) +BHN*GNM2T (I) ) 
1  /SM(I) 

1180  CONTINUE 
C  TRANSFORM  DELT  TO  DEL 
1=0 

DO  1220  IJ=1,NUP 
KJ=K+(IJ-1 )*NUMNP 
DO  1220  L=1 ,2 

II=(K-1 )*i+L+(IJ-1 )*NUMNP*3 

DEL(II)=0.0 

1=1+1 

ID=CODE(KJ ) 

IF(ID.EQ.O)GO  TO  1195 

IF( (ID. EQ. 1 .OR. ID. EQ. 6) .AND.L.EQ. 1 )G0  TO  1185 
IF( (ID.EQ.2.0R.ID.EQ.5) .AND.L.EQ.2)G0  TO  1185 
IF ( ID . EQ . 4 . OR . ID . EQ . 7 ) GO  TO  1185 
IF( (ID.EQ.8.0R.ID.EQ. 13) .AND.L.EQ. 1 )G0  TO  1190 
IF ( ( ID . EQ . 9 . OR . ID . EQ . 1 2 ) . AND . L . EQ . 2 )G0  TO  1190 
IFdD.EQ.11  .0R.ID.EQ.14)G0  TO  1190 
GO  TO  1195 
1185  CALL  BCDEL(II) 

GO  TO  1220 
1190  CALL  BCVEL(II) 

GO  TO  1220 
1195  DO  1200  J=1 ,4 

1200  DEL(II )=DEL(II)+TF(I , J)*DELT( J ) 

1220  CONTINUE 
RETURN 

1230  WRITE (6 , 1250) 

1250  FORMAT ("  *  *  *  »  «  ERROR  IN  SUBROUTINE  BAND  **•••«/ 

1  "  *  *  *  *  *  CHECK  KO  OR  BOUNDARY  CONDITIONS  »*»»»") 

END 

SUBROUTINE  INTER 
INTEGER  CODE 


87 


COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/MATP/RO( 12) ,E(9 , 12) ,EE(9 ) ,ETA( 12) 

COMMON /ARG /XXX ( 10 ) , YYY (10),S(24) ,XX(3),YY(3) , 

1 CRZ (6, 6), XI (IQ), SIG ( 1 2 ) ,  N ,  M 

COMMON/ELDATA/BETA (12), ALPHA ( 1 2 ) , TH ( 1 2 ) , IX ( 200 , 4 ) , MATRIL (12) 

COMMON /RESULT/D ( 6, 6 ),C( 6, 6) ,CNS(6,6) 

COMMON/TD/IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAXJ , NMTL , NBC 
COMMON/BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NT1ME 
COMMON /SIGM/BSUM , TSUM ( 6 ) , SIGMA ( 1 2 , 50 , 6 ) , DSIG ( 6 ) , F ( 60  0 ) 

COMMON /DISP/DELN 1 (600 ) ,DELN(600 ) ,DEL(600 ) ,GNM1 (600 ) ,GNM2(600 ) 
COMMON/MASS/A(600) ,B(600) ,CM(8) 

DIMENSION  X(7),Y(7),QQ(9) 

DATA  QQ/3* . 1259391 805448 , 3* . 1 323941527884 , .225 , 

1  .696140478028, .410426192314/ 

X ( 7 ) = (XX ( 1 )+XX(2 )+XX(3) )/3 .0 
Y(7)=(YY(1 )+YY(2)+YY(3))/3.0 
DO  100  1=1,3 
J=I+3 

X(I)=QQ(8)»XX(I)+(1.00-QQ(8))«X(7) 

X(J)=QQ(9) *XX ( I ) + ( 1 . OO-QQ ( 9 ) ) *X ( 7 ) 

Y(I) =QQ(8  )*YY (I )  +  ( 1 .00-QQ(8 ) )*Y(7 ) 

100  Y(J)=QQ(9)*YY(I)+(1 ,00-QQ(9))*Y(7 ) 

DO  300  1=1,10 
300  XI(I)=0.00 

AREA=.50«(XX(1 )»(YY(2)-YY(3))+XX(2)»(YY(3)-YY(1 ))+XX(3)*(YY(1 ) 

1  -YY ( 2 ) ) ) 

DO  400  1=1,7 
XI  (1  )=XI(1 )+QQ(I) 

XI(2)=XI(2) +QQ ( I ) *X ( I ) 

XI ( 3 ) =Xl ( 3 ) +QQ ( I ) *Y ( I ) 

XI(4)=XI(4)  +QQ (I)*X(I)*Y(I) 

XI(5)=XI(5)+QQ(I)»X(I)»»2 
400  XI (6 )=XI(6 )+QQ(I )*Y(I )**2 
DO  500  1=1,10 
500  XI (I )=XI(I)*AREA 
RETURN 
END 

SUBROUTINE  LOAD 
INTEGER  CODE 

COMMON /BASIC/ VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/ARG/XXX( 10 ) ,YYY ( 10 ) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 ,6) ,XI(10) ,SIG(12) ,N,M 
COMMON /NPDATA/X ( 200 ) , Y ( 200 ) , CODE ( 200 ) , NPNUM( 1 0 , 20 ) 

COMMON /ELDATA/BETA (12) , ALPHA ( 12),TH(12) , IX (200 , 4 ) ,MATRIL( 1 2) 
COMMON /BASIC2 /BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM , TSUM( 6 ) , SIGMA( 12,50,6), DSIG (6) ,F(600) 

COMMON /FOR /FF (5) ,FC(5) 

DO  50  1  =  1, IT 
F(I )=0 .00 
50  CONTINUE 
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C  NUMBER  NODES  OF  THE  RECTANGULAR  ELEMENTS 

DO  20  N=1 ,NUMEL 
I1=IX(N, 1 ) 

I2=IX(N ,2) 

I3=IX(N,3) 

I4=IX(N,4) 

C  DESIGNATE  COORDINATES 

XXX(1)=X(I1) 

XXX(2)=X(I2) 

XXX(3)=X(I3) 

XXX(4)=X(I4) 

XXX(5 )=1 .  0/4 .0*(XXX( 1 )+XXX(2)+XXX (3)+XXX(4 ) ) 

YYY (1)=Y(I1) 

YYY(2 )=Y(I2) 

YYY(3)=Y(I3) 

YYY(4)=Y(I4) 

YYY(5 )=1 .0/4 .0*( YYY ( 1 )+YYY(2)+YYY(3 )+YYY(4 ) ) 

CALL  DISFOR 

IF ( FF ( 1 ) . EQ . 0 . 0 . AND . FF ( 2 ) . EQ . 0 . 0 . AND . FF ( 3 ) • EQ . 0 . 0 . AND . FF ( 4 ) . EQ . 0 . 0 
1)  GO  TO  20 

C  CALL  LOADS  OF  TRIANGLES 

DO  9  1=1,5 
FC(I)=0.00 
9  CONTINUE 
CALL  LOT( 1,2) 

CALL  LOT (2, 3) 

CALL  LOT (3 ,4) 

CALL  LOT (4,1) 

DO  21  1=1,4 
FC(I)=FC(I)+FC(5)/4.0 
21  CONTINUE 

C  CHANGE  TO  GLOBAL  FORCES  IN  W-DIRECTION 

II 1 =NUMLA*NUMNP*3+i 1 *3 
II2=NUMLA*NUMNP*3+I2* 3 
II3=NUMLA*NUMNP*3+I3*3 
II4=NUMLA*NUMNP *3+14*3 
F(II1 )=F(II1 )+FC( 1 ) 

F(II2)=F(II2)+FC(2) 

F(II3)=F(II3)+FC(3) 

F(II4)=F(II4)+FC(4 ) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  LOT(II,JJ) 

INTEGER  CODE 

COMMON/BASIC/ VOL , NUMNP , NUMEL , NUMLA , NCG 
C0MM0N/ARG/XXX(10),YYY(10) ,S(24),XX(3) ,YY(3) , 

1CRZ(6 , 6) ,XI( 10) ,SIG( 12)  ,N,M 

COMMON /NPDATA/X ( 200 ) , Y ( 200 ) , CODE ( 200 ) , NPNUM ( 1 0 , 20 ) 

COMMON /ELDAT A/BETA ( 1 2 ) , ALPHA ( 1 2 ) , TH ( 1 2 ) , IX ( 200 , 4 ) , MATRIL ( 1 2 ) 
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COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM , TSUM( 6 ) , SIGMA (12,50,6),DSIG(6),F(600) 

COMMON /FOR /FF(5 ) ,FC(5) 

DIMENSION  DD1 (3) ,DD2(3) ,DD3(3) ,DD(3,3) ,FE(3) ,FFF(j) 

DIMENSION  AA(3)»BB(3) »CC(3) 

C  DESIGNATE  THE  TRIANGULAR  DISTRIBUTED  FORCES 
FFF(1)=FF(II) 

FFF(2)=FF(JJ) 

FFF(3)=FF(5) 

XX(1)=XXX(II) 

XX(2)=XXX(JJ) 

XX(3)=XXX(5) 

YY( 1 )=YYY (II) 

YY (2 )=YYY ( JJ ) 

YY (3)=YYY (5 ) 

AA(1 )=XX(2)»YY(3)-XX(3)*YY(2) 

AA(2 ) =XX (3 )*YY ( 1 )-XX(1 )*YY(3) 

AA(3)=XX(1)»YY(2)-XX(2)»YY(1) 

BB( 1 )=YY (2)-YY (3 ) 

BB(2)=YY(3)-YY(1) 

BB(3)=YY(1 )-YY(2) 

CC ( 3 ) =XX ( 2 ) -XX  ( 1 ) 

CC ( 2 ) =XX ( 1 ) -XX ( 3 ) 

CC  ( 1 )  =XX  ( 3 )  -XX  ( 2 ) 

C  INTEGRATE  XX  AND  YY 

CALL  INTER 
DO  12  1=1 ,3 

DD1 (I)=AA(I)*XI(1 )+BB(I)*XI(2)+CC(I)*XI(3) 

DD2(I)=AA(I)*XI(2) +BB (I)*XI(5) +CC (I)*XI(4) 
DD3(I)=AA(I)*XI(3)+BB(I)«XI(4)+CC(I)»XI(6) 

12  CONTINUE 
DO  18  1=1,3 
DO  18  J=1 ,3 

DD(I,J)=AA(I)*DD1(J)+BB(I)*DD2(J)+CC(I)*DD3(J) 

18  CONTINUE 

C  CALCULATE  EQUIVELENT  CONCENTRATED  FORCES 

AREA=.50*(XX(1)*(YY(2)-YY(3))+XX(2)»(YY(3)-YY(1))+XX(3)*(YY(1) 

1  -YY(2 ) ) ) 

DO  99  1=1,3 

FE(I )= 1 .0/(4 .0*AREA**2)*(DD(I , 1 )*FFF( 1 )+DD(I , 2)*FFF(2 )+DD(I , 3)*FFF 
1(3)) 

99  CONTINUE 

FC(II)=FC(II)+FE(  1 ) 

FC  ( J  J )  =FC  ( J  J )  +FE  ( 2 ) 

FC(5 )=FC(5 )+FE(3) 

RETURN 

END 

SUBROUTINE  MESH 
INTEGER  CODE 

DIMENSION  AR( 1 0 ,  40),AZ(10,  40)  ,NC0DE(10,  40) 
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COMMON /BASIC/VOL ,NUMNP ,NUMEL,NUMLA ,NCG 
COMMON /MAT  P /RO ( 1 2 ) ,E(9 , 12) ,EE(9) ,ETA(12) 

COMMON / ARG /XXX (10)  ,m(10)  ,S(24)  ,XX(3)  ,YY(3) , 

1CRZ(6,6) ,XI(10) ,SIG(12) ,N,M 
COMMON /NPDATA/X (200 ) ,Y(200 ) , CODE (200 ) ,NPNUM( 10 , 20) 
C0MM0N/ELDATA/BETA(12) , ALPHA (12) ,TH(12) ,IX(200 ,4) ,MATRIL(12) 
COMMON/ RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

COMMON /TD /IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
C0MM0N/BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM,TSUM( 6) ,SIGMA( 12 , 50 , 6) ,DSIG(6) ,F(600) 
C0MM0N/DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/ A(600) ,B(600) ,CM(8) 

EQUIVALENCE  (X(1 ) ,AR) , (Y (1 ) ,AZ) ,(IX(1 , 1 ) ,NCODE) 

QXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  MESH  CONTROL  INFORMATION 

QXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

READ  (5,1000)  MAXI, MAX J,NSEG, NBC, NMTL 

WRITE (6 , 2000 )  MAXI , MAX J , N3EG , NBC , NMTL 

QX  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

C  INITIALIZE 

QXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

ISEG=-1 
PI=3. 1415927 
DO  110  J=1 ,40 
DO  100  1=1,10 
NC0DE(I,J)=0 
AR(I,J)=0. 

AZ(I , J )=0 . 

JMAX(I )=0 
100  JMIN(I)=MAXI 
IMIN ( J )=MAXJ 
110  MX(J)=0 

QXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

C  LINE  SEGMENT  CARDS 

c*  *********************************** 

150  ISEG=ISEG+1 

159  IF ( ISEG . EQ . NSEG )  GO  TO  400 

READ (5, 1001 )  II ,J1 ,R1 ,Z1 ,12,J2,R2,Z2,I3,J3,R3,Z3,IPTION 

WRITE (6 ,2001 )I1 ,J1 ,R1 ,Z1 ,I2,J2,R2,Z2,I3,J3,R3,Z3,IPTI0N 

IPTI0N=IPTI0N+1 

AR(I1 ,J1 )=R1 

AZ(I1,J1)=Z1 

NC0DE(I1 , J 1 ) = 1 

CALL  MNIMX(I 1 , J 1 ) 

GO  TO  (150,200,200,300,300,200,200),  IPTION 

C*  *********************************** 

C  GENERATE  STRAIGHT  LINES  ON  BOUNDARY 

Q*  *********************************** 

200  DI=  ABS(FLOAT (12-11 )) 

DJ=  ABS(FLOAT (J2-J1 ) ) 
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AR(I2,J2)=R2 
AZ(I2,J2)=Z2 
NC0DE(I2,J2)=1 
CALL  MNIMX(I2,J2) 

ISTRT=I1 

ISTP=I2 

JSTRT=J1 

JSTP=J2 

DIFF=MAX1 (DI,DJ) 

ITER=DIFF-1 . 

IINC=0 

JINC=0 

IF(I2.NE.I1)  IINC=(I2-I1 )/IABS(I2-I1 ) 

IFCJ2.NE.J1 )  JINC=(J2-J1 )/IABS(J2-J1 ) 

KAPPA=1 

IF(I2.NE.I1 .AND. J2.NE.J1 .AND.IPTI0N.NE.3)  KAPPA=2 
IFCKAPPA.EQ.2)  DIFF=2.»DIFF 
RINC= (R2-R1 )/DIFF 
ZINC=(Z2-Z1 )/DIFF 

WRITE ( 6 , 2002 )  DI ,  DJ , DIFF , RINC , ZINC , ITER , IINC , JINC , KAPPA 

CHECK  FOR  INPUT  ERROR 

IF ( KAPPA. NE. 2. OR .DI .EQ.DJ)  GO  TO  210 
WRITEC6 , 2003 ) 

GO  TO  150 

INTERPOLATE 


210  1=11 
J=J1 

WRITEC6 ,2004) 

DO  230  M=1 , ITER 

IFCITER.EQ.0.AND.IPTI0N.EQ.2)  GO  TO  230 

IF ( ITER . EQ . 0 . AND . IPTION . EQ . 6 )  GO  TO  230 

IF ( ITER . EQ . 0 . AND . IPTION . EQ . 7 )  GO  TO  230 

IFCKAPPA.EQ.2)  GO  TO  220 

IOLD=I 

I=I+IINC 

JOLD=J 

J=J+JINC 

AR ( I , J ) = AR ( IOLD , JOLD ) +RINC 
AZ ( I , J ) =AZ ( IOLD ,  JOLD ) +ZINC 
WRITEC6 ,2005)  I,J,AR(I,J) ,AZ(I, J) 

CALL  MNIMXCl , J) 

NCODECI, J)=1 
GO  TO  230 
220  CONTINUE 

IF(I1 .GT . 12. AND. IPTION .EQ. 7 )  GO  TO  221 
IFCII.LT. 12. AND. IPTION. EQ. 6)  GO  TO  221 
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IOLD=I 

I=I+IINC 

AR (I , J ) =AR (IOLD , J )+RINC 
AZ ( I , J ) =AZ ( IOLD , J )+ZINC 
WRITE(6,2005J  I ,  J ,AR(I , J) ,AZ(I , J) 

N CODE (I , J)=1 
CALL  MNIMX(I , J) 

JOLD=J 

J=J+JINC 

AR ( I , J ) =AR ( I , JOLD )+RlNC 
AZ ( I , J ) = AZ ( I , JOLD )+ZINC 
NCODE(I , J )= 1 

WRITE(6 ,2005)  I , J , AR (I , J ) , AZ ( I , J ) 

CALL  MNIMX(I , J) 

GO  TO  230 
221  JOLD=J 
J=J+JINC 

AR(I , J )=AR(I , JOLD)+RINC 
AZ(I,J)=AZ(I , JOLD)+ZINC 
NCODE(I , J)=1 

WRITE(6 ,2005)  I,J,AR(I,J),AZ(I,J) 

CALL  MNIMX(I , J) 

IOLD=I 

I=I+IINC 

AR (I , J ) = AR (IOLD , J )+RINC 
AZ(I,J)=AZ(IOLD,J)+ZINC 
NCODE(I , J )=1 

WRITE(6 ,2005)  I, J,AR(I,J) ,AZ(i,J) 

CALL  MNIMX(I,J) 

230  CONTINUE 

IF(KAPPA.EQ.I)  GO  TO  150 
IF(I1 .GT.I2.AND.IPTI0N.EQ.7)  GO  TO  231 
IF ( I 1 . LT . 12 . AND . IPTION . EQ . 6 )  GO  TO  231 
IOLD=I 
I=I+IINC 

AR (I , J ) =AR ( IOLD , J )+RINC 
AZ(I , J )=AZ(IOLD , J)+ZINC 
GO  TO  232 

231  CONTINUE 
JOLD=J 
J=J+JINC 

AR(I,J)=AR(I,JOLD)+RINC 

AZ(I,J)=AZ(I,JOLD)+ZINC 

232  CONTINUE 
NCODE(I ,  J  )  =  1 

irfRITE(6 ,2005)  I,J,AR(I,J),AZ(I,J) 

CALL  MNIMX(I , J) 

GO  TO  150 

c*  *********************************** 

C  GENERATE  CIRCULAR  ARCS  ON  BOUNDARY 
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DI=  ABS ( FLOAT ( ISTP-ISTRT ) ) 

DJ=  ABS ( FLOAT ( JSTP-JSTRT) ) 

IINC=0 

JINC=0 

IF ( ISTRT . NE . ISTP )  IINC=  ( ISTP-ISTRT ) /IABS ( ISTP-ISTRT ) 
IF ( JSTRT . NE . JSTP )  JINC= ( JSTP-JSTRT ) /IABS ( JSTP-JSTRT ) 
LAMDA=1 

IF(IINC.NE.O.AND. JINC.NE.O)  LAMDA=2 
DIFF=MAX1(DI,DJ) 

ITER=DIFF-1 . 

IF(LAMDA.EQ.2)  DIFF=2.«DIFF 

DELPHI=(ANG2-ANG1 )/DIFF 

WRITE (6, 2008)  ANG1 ,ANG2,DIFF, DELPHI 

CHECK  FOR  INPUT  ERROR 

IF ( LAMDA . NE . 2 . OR . DI . EQ . DJ )  GO  TO  350 
WRITE(o , 2003 ) 

GO  TO  150 
350  IO=ISTRT 
JO=JSTRT 
WRITE (6, 2004) 

INTERPOLATE 

NPT=IABS (12-11 )+IABS(J2-J1 )-1 
DO  380  M=1,ITER 

359  IF(LAMDA.EQ.2)  GO  TO  360 
I=IO+IINC 

J=JO+JINC 
CALL  MNIMX(I, J) 

NCODE(I, J)=1 

CALL  CIRCLE ( ANG 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 

WRITE (6, 2005)  I,J ,AR(I,J) ,AZ(I,J) 

GO  TO  370 

360  I=IO+IINC 
J=JO 

NCODE(I , J)= 1 
CALL  MNIMX(I,J) 

CALL  CIRCLE ( ANG 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 
WRITE(6 ,2005 )  I,J,AR(I,J),AZ(I,J) 

J=JO+JINC 
NC0DE(I,J)=1 
CALL  M1'JIMX(I, J) 

CALL  CIRCLE ( ANG 1 , DELPHI , RSTRT , ZSTRT , RC , ZC , I , J ) 
WRITE(6 ,2005)  I,J,AR(I,J) ,AZ(I,J) 

370  10=1 
380  JO=J 

IF(LAMDA.NE.2)  GO  TO  390 
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C*  ****«******************««»«*«««**** 

300  Afi(I2,J2)=R2 
AZ(I2, J2)=Z2 
NCODE(I2,J2)  =  1 
CALL  MNIMX(I2 ,'J2 ) 

IF ( IPTION . EQ . 5 )  GO  TO  320 

FIND  CENTER  OF  CIRCLE 

AR(I3,J3)=R3 
AZ(I3,J3)=Z3 
NCODE(I3,J3)=1 
CALL  MNIMX(I3,J3) 

SLAC=(Z2-Z1 )/(R2-R1 ) 

SLBF=-1 ./SLAC 
SLCE=(Z3-Z2)/(R3-R2) 

SLDF=-1 . /SLCE 

CHECK  FOR  INPUT  ERROR 

IF ( ABS ( SL AC-SLCE ) . GT . .001 )  GO  TO  310 
WRITE(6 ,2006 )  R1 ,Z1 ,R2,Z2,R3,Z3, SLAC, SLCE 
GO  TO  150 

310  R4=R1+(R2-R1 )/2. 

Z4=Z1+(Z2-Z1 )/2 . 

R5=R2+(R3-R2)/2. 

Z5=Z2+(Z3-Z2)/2. 

BBF=Z4-SLBF*R4 
BDF=Z5-SLDF*R5 
RC= ( BBF-BDF ) / ( SLDF-SLBF ) 

ZC=SLBF*RC+BBF 
WRITE(6 ,2007 )  RC , ZC 
KAPPA=1 
GO  TO  330 
320  KAPPA=2 
RC=R3 
ZC=Z3 

330  ISTRT=I1 
ISTP=I2 
JSTRT=J1 
J3TP=J2 
RSTRT=R1 
RSTP=R2 
Z3TRT=Z1 
ZSTP=Z2 

340  CALL  ANGLE ( R3TRT, ZSTRT, RC, ZC, ANG 1 ) 

CALL  ANGLE ( RSTP , ZSTP , RC , ZC , ANG2 ) 

IF(ANG2.LE.ANG1 )  ANG2=2 . 0*PI+ANG2 
C 

C  FIND  ANGULAR  INCREMENT 
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I=IO+IINC 
NCODE(I , J)=1 
CALL  MNIMX ( I , J ) 

CALL  CIRCLE ( ANG 1 , DELPHI , R3TRT , ZSTRT ,RC,ZC,I, J) 

WRITE (6 ,2005)  I,J,AR(I, J) ,AZ(I,J) 

390  IFCKAPPA.EQ.2)  GO  TO  150 
ISTRT=I2 
ISTP=I3 
JSTRT=J2 
JSTP=J3 
RSTRT=R2 
RSTP=R3 
ZSTRT=Z2 
ZSTP=Z3 
KAPPA=2 

399  GO  TO  340 

c»  a********************************** 

C  CALCULATE  COORDINATES  OF  INTERIOR  POINTS 

C*  *********************************** 

400  IF(MAXJ.LE.2)  GO  TO  430 
J2=MAXJ-1 

DO  420  N= 1,500 
RESID=0 . 

DO  410  J=2,J2 
I 1 =IMIN ( J)+1 
I2=IMAX(J)-1 
DO  410  1=11,12 

IF(NC0DE(I,J).EQ.1)  GO  TO  410 

DR=(AR(I+1 ,J)+AR(I-1 , J )+AR(I , J+1 )+AR(I , J-1 ) ) /4 . -AR(I , J) 

DZ=(AZ(I+1 , J)+AZ(I-1 ,J)+AZ(I,J+1)+AZ(I,J-1))/4.-AZ(I,J) 
RESID=RESID+AB3(DR)+ABS(DZ) 

AR(I,J)=AR(I ,J)+1 ,8*DR 
AZ(I,J)=AZ(I,J)+1.8*DZ 
410  CONTINUE 

IF(N.EQ.I)  RES1 =RESID 
IF(N .EQ. 1 .AND.RESID.EQ.O. )G0  TO  430 
IF(RESID/RES1.LT.1.E-5)  GO  TO  430 
420  CONTINUE 
430  WRITE(6 ,2009)  N 

Q*  *********************************** 

CALL  POINTS 

c*  *********************************** 

1000  FORMAT  (515) 

1001  FORMAT  (3(2X3, 2F 8.3) ,15) 

2000  FORMAT  (30H1  MESH  GENERATION  INFORMATION// 


1  41  HO  MAXIMUM  VALUE  OF  I  IN  THE  MESH - 13/ 

2  41  HO  MAXIMUM  VALUE  OF  J  IN  THE  MESH - 13/ 

3  41H0  NUMBER  OF  LINE  SEGMENT  CARDS - 13/ 

4  41 HO  NUMBER  OF  BOUNDARY  CONDITION  CARDS - 13/ 

5  41  HO  NUMBER  OF  MATERIAL  BLOCK  CARDS - 13///) 
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2001  FORMAT  (//88H  INPUT  II  J1  R1  Z1  12  J2  R2  Z 

12  13  J3  R3  Z3  IPTION/8x,3(2I4,2F8.4) ,16) 

2002  FORMAT  (5H  DI=F4.0,5H  DJ=F4.0,7H  DIFF=F4.0,7H  RINC=F8.3,7H  ZI 

1NC=F8.3,7H  ITEfl=I3,7H  IINC=I3,7H  JINC=I3,8H  KAPPA=I1) 

2003  FORMAT  ( IX,  38H»*BAD  INPUT— THIS  LINE  IS  NOT  DIAGONAL) 

2004  FORMAT  (30H  I  J  AR  AZ) 

2005  FORMAT  (2I5.2F11.6) 

2006  FORMAT  (51 H  **  BAD  INPUT  -  THESE  POINTS  DO  NOT  DEFINE  A  CIRCLE,/, 
13X, 6F12 .4 , 10X , 2E20 .8) 

2007  FORMAT ( 1 9H  CENTER  COORDINATE, (F1 1 .6 , 1X,F1 1 .6 , IX) ) 

2008  FORMAT  (7H  ANG1=F9.6,7H  ANG2=F9-6,7H  DIFF=F3.0,9H  DELPHI =F9. 6) 

2009  FORMAT  (//30H  COORDINATES  CALCULATED  AFTER  13,11H  ITERATIONS) 
RETURN 

END 

SUBROUTINE  MNIMX(I,J) 

INTEGER  CODE 

COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /MAT P/RO( 12) ,E(9, 12) ,EE(9) ,ETA(12) 

COMMON / ARG/XXX (10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 ,6) ,XI( 10) ,SIG(12) ,N ,M 

COMMON /NPDATA/X (200 ) ,Y(200) .CODE (200) ,NPNUM( 10 ,20) 

COMMON/ELD ATA/BETA (12), ALPHA ( 1 2 ) , TH ( 1 2 ) , IX ( 200 , 4 ) , MATRIL (12) 

COMMON /RESULT /D( 6, 6) ,C(6 ,6) ,CNS(6 , 6) 

COMMON /TD /IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
COMMON /BASI C2 /BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NL AY , IFLAG , IT , NTlME 
COMMON /SIGM/BSUM , TSUM ( 6 ) , SIGMA (12,50,6) ,DSIG(6) ,F(600) 

COMMON /DISP/DELN 1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/A(600) ,B(600) ,CM(8) 

IF ( J . LT . JMIN ( I ) )  JMINd)=J 
IF ( J . GT . JMAX ( I ) )  JMAX(I )=J 
IF ( I . LT . IMIN ( J ) )  IMIN(J)=I 
IF(I.GT.IMAX(J))  IMAX( J )=I 
RETURN 
END 

SUBROUTINE  MPLOT 
INTEGER  CODE 

COMMON /BASI C/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /MATP/RO( 12 ),E(9, 12) ,EE(9) ,ETA( 12) 

COMMON /ARG /XXX ( 10 ) ,YYY (10) ,S(24)  ,XX(3)  ,YY(3), 

1CRZ(6 ,6) , XI ( 1 0 ) , SIG (12) , N , M 
COMMON /NPD ATA/R ( 200 ) , Z ( 200 ) , CODE ( 200 ) , NPN  UM ( 1 0 , 20 ) 

COMMON /ELDAT A/BETA (12), ALPHA (12) ,TH( 12) , IX (200 , 4) ,MATRIL(12) 

COMMON /RESULT/D (6, 6) ,C(6, 6) ,CNS(6, 6) 

COMMON /TD /IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
C0MM0N/BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON/SIGM/BSUM , TSUM( 6 ) , SIGMA (12,50,6), D3IG ( 6 ) , F ( 600 ) 
COMMON/DISP/DELN 1(600) ,DELN(600) , DEL (600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/A (600) ,B(600) ,CM(8) 

DIMENSION  TITLE (8) ,THI( 13) 

IF (NTIME. GT.O)GO  TO  50 
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READ (5,100) TITLE , RMAX , ZMAX 
C  WRITE (6, 100) TITLE, RMAX, ZMAX 
TITLE (8)="; " 

50  CONTINUE 

100  FORMAT (8A10/2F1 0.5) 

NUP=NUMLA+1 

S3=.5 

C3=SQRT(3.)/2. 

THI(1 )=0.0 
DO  200  1=1  ,NUMLA 
200  THI(I+1 )=TH(I)+THI(I) 

XMAX=C3*RMAX»1 .1 
XMIN=-C3*ZMAX*1 .1 
YMAX= (XMAX-XMIN )*3 . /4 . 

YMIN=-YMAX/3. 

C  DRAW  X  LINES 

CALL  UWIN DO ( XMIN , XMAX , YMIN , YMAX ) 
YMIN=YMIN/4 . 

CALL  UPRINT(XMIN,0. .TITLE) 

CALL  UPRINT (XMIN, YMIN,"  TIME  =  ;") 
CALL  UPRNT1 (TIME, "REAL") 

DO  300  J=1 ,MAXJ 
NSTART=IMIN(J) 

NSTOP=IMAX( J ) 

DO  300  K=1 ,NUP 

IF(J.GT. 1 .AND.K.NE.NUP)GO  TO  300 
DO  250  I=NSTART,NST0P 
NP=NPNUM(I,J) 

X=C3*(R(NP)-Z(NP) ) 
Y=THI(K)+S3*(R(NP)+Z(NP) ) 

IF ( I . EQ . NSTART ) CALL  UMOVE(X,Y) 

CALL  UPEN(X,Y) 

250  CONTINUE 
300  CONTINUE 
C  DRAW  Y  LINES 

DO  400  1=1, MAXI 
NSTART=JMIN (I ) 

NSTOP=JMAX(I ) 

DO  400  K=1 ,NUP 

IF(I ,GT . 1 .AND .K.NE.NUP )G0  TO  400 
DO  350  J  =NSTART , NSTOP 
NP=NPNUM(I , J) 

X=C3*(R(NP)-Z(NP) ) 

Y=THI (K)+S3*(R(NP )+Z(NP ) ) 

IF (J.EQ. NSTART) CALL  UMOVE(X,Y) 

CALL  UPEN(X,Y) 

350  CONTINUE 
400  CONTINUE 
C  DRAW  Z  LINES 

DO  500  1=1, MAXI 
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NP=NPNUM(I , 1 ) 

DO  500  K=1 ,NUP 
X=C3*(fi(NP)-Z(NP) ) 

Y=S3*(R(NP)+Z(NP))+THI(K) 

IF(K.EQ. 1 )CALL  UMOVE(X,Y) 

500  CALL  UPEN (X , Y ) 

DO  600  J=1,MAXJ 
NP=NPNUM( 1 , J ) 

DO  600  K=1 ,NUP 
X=C3*(R(NP)-Z(NP) ) 

YrS3*(R(NP)+Z(NP))+THI(K) 

IF(K.EQ. 1 )CALL  UMOVE(X,Y) 

600  CALL  UPEN (X , Y) 

IF(NTIME.EQ.O)GO  TO  650 
CALL  MPLOTD ( ZMAX , RMAX ) 

GO  TO  800 
C  LABEL  NODES 

650  CALL  USET( "INTEGER") 

DO  700  1=1, MAXI 
DO  700  J  =  1  ,  MAX  J 
DO  700  K=1 , NUP 
NP=NPNUM(I, J) 

X=C3*(R(NP)-Z(NP) ) 

Y=S3*(R(NP)+Z(NP))+THI(K) 

PI =NP+(K-1 )*NUMNP 

700  IF ( I . EQ . 1 . OR . J . EQ . 1 . OR . K . EQ . NUP ) CALL  UPRINT(X,Y,P1 ) 

CALL  USETC'TEXT1') 

800  CALL  UERASE 
RETURN 
END 

SUBROUTINE  MPLOTD ( ZMAX, RMaX) 

INTEGER  CODE 

COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /MATP/RO(1 2) ,E(9 ,12) ,EE(9) ,ETA(12) 

COMMON /ARG/XXX( 10) ,YYY(10) ,S(24) ,XX(3),YY(3) , 

1CRZ(6 ,6) ,XI(10) ,SIG(12)  ,N,M 
COMMON/NPDATA/R(200 ) ,Z(200 ) ,CODE(200) ,NPNUM( 10 , 20 ) 

COMMON /ELDATA/BETA( 12) , ALPHA ( 12) ,TH(12) ,IX(200,4) ,MATHIL(12) 

COMMON /RESULT/D (6, 6) ,C(6,6) ,CNS(6,6) 

COMMON/TD/IMIN ( 100 ) ,IMAX( 100 ) , JMIN (25 ) , JMAX(25 ) , MAXI , MAX J ,NMTL ,NBC 
COMMON /BA3IC2/ BET ,H,HN,HH(3),HT(2), TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM, TSUM(6 ), SIGMA (12, 50, 6) ,DSIG (6 ) ,F(600 ) 

COMMON /DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600 ) ,GNM2(600 ) 
COMMON /MASS/ A (600 ) ,B(600 ) ,CM(8 ) 

DIMENSION  TITLE (8) ,THI(13) 

CALL  USETC'DASH") 

NUP=NUMLA+1 

S3=.5 

C3=SQRT(3.)/2. 

THI ( 1 ) =0 .0 
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DO  200  1=1 ,NUMLA 
200  THI(I+1 )=TH(I)+THI(I) 

C  DRAW  X  LINES 

DO  300  J = 1 , MAX J 
NSTART=IMIN(J) ' 

NSTOP=IMAX( J) 

DO  300  K= 1 ,NUP 

IF( J.GT. 1 ,AND.K.NE.NUP)GO  TO  300 
DO  250  I=NSTART,NSTOP 
NP=NPNUM(I , J) 

NP1=(NP+(K-1 )*NUMNP)*3 
X=C3*(R(NP)+DEL(NP1-2)-Z(NP)+DEL(NP1-1)) 

Y=THI(K)+DEL(NP1 )+S3*(R(NP)+DEL(NP1-2)+Z(NP)+DEL(NP1-1  ) ) 
IF ( I . EQ . NSTART ) CALL  UMOVE(X,Y) 

CALL  UPEN(X,Y) 

250  CONTINUE 
300  CONTINUE 
C  DRAW  Y  LINES 

DO  400  1=1, MAXI 
NSTART= JMIN (I ) 

NSTOP=JMAX(I ) 

DO  400  K=1 ,NUP 

IF(I.GT.1 .AND.K.NE.NUP)GO  TO  400 
DO  350  J=NSTART,NSTOP 
NP=NPNUM(I , J) 

NP1 = (NP+(K-1 )*NUMNP)*3 
X=C3*(R(NP)+DEL(NP1-2)-Z(NP)+DEL(NP1-1 )) 
Y=THI(K)+DEL(NP1)+S3*(R(NP)+DEL(NP1-2)+Z(NP)+DEL(NP1-1 )) 
IF (J.EQ. NSTART) CALL  UMOVE(X,Y) 

CALL  UPEN(X,Y) 

350  CONTINUE 
400  CONTINUE 
C  DRAW  Z  LINES 

DO  500  1=1, MAXI 
NP=NPNUM(I , 1 ) 

DO  500  K=1 ,NUP 
NP1=(NP+(K-1 )*NUMNP)*3 
X=C3*(R(NP)+DEL(NP1-2)-Z(NP)+DEL(NP1-1)) 
Y=S3*(R(NP)+DEL(NP1-2)+Z(NP)+DEL(NP1-1 ))+THI(K)+DEL(NP1) 
IF(K.EQ. 1 )CALL  UMOVE(X,Y) 

500  CALL  UPEN (X , Y) 

DO  600  J=1 ,MAXJ 
NP=NPNUM( 1 , J ) 

DO  600  K=1 ,NUP 
NP1=(NP+(K-1 )*NUMNP )*3 

X=C3*(R(NP)+DEL(NP1-2)-Z(NP)+DEL(NP1-1 )) 
Y=S3*(R(NP)+DEL(NP1-2)+Z(NP)+DEL(NP1-1 ))+THI(K)+DEL(NP1) 
IF(K.EQ.1 )CALL  UMOVE(X,Y) 

600  CALL  UPEN(X,Y) 

CALL  USET( "LINE") 
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RETURN 

END 

FUNCTION  NODE(I ,  J ) 

INTEGER  CODE 

COMMON  /BASIC/VO'L ,  NUMNP ,  NUMEL ,  NUMLA ,  NCG 
COMMON /MAT P/RO( 12) ,E(9 , 12 ) ,EE(9 ) ,ETA( 12) 

COMMON /ARG /XXX (10)  ,YYY  (10)  ,S(24)  , XX ( 3 )  ,W(3)  , 

1CRZ(6 ,6) , XI (10)  ,SIG(12) ,N,M 
C0MM0N/NPDATA/X(200 ) , Y (200 ) , CODE (200 ) ,NPNUM( 1 0 , 20 ) 
COMMON/ELDATA/BETA( 12) ,ALPHA(12) ,TH( 12) ,IX(200 ,4) ,MATRIL(12) 
COMMON/RESULT/D(6,6) ,C(6,6) ,CNS(6,6) 

COMMON /TD/IMIN ( 1 00 ) ,IMAX( 100) , JMIN (25 ) , JMAX(25 ) , MAXI , MAX J ,NMTL , NBC 
COMMON/BASIC2/BET ,H,HN,HH(3) ,HT(2) , TIME ,NLAY , IFLAG , IT ,NTIME 
COMMON/SIGM/BSUM,TSUM(6) ,SIGMA(12,50,6) ,D3IG(6) ,F(600) 

COMMON /DISP/DELN 1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON/MASB/A(600) ,B(600) ,CM(8) 

NODfi=0 

DO  100  JJ=1  ,J 
NSTART=IMIN ( J J ) 

NSTOP=IMAX( JJ ) 

DO  100  II =NSTART , NSTOP 
NODE=NODE+1 

IF(JJ.EQ.J. AND . II . EQ . I )  RETURN 
100  CONTINUE 
RETURN 
END 

SUBROUTINE  POINT (TIME , NUM , TN ) 

COMMON /NUMB /NITER , HD AT , TO 
AN= (TIME-TO) /HDAT+1 .0 
NUM=AN 

NUM= ( AN+NUM+1 .5)/2 
TN=TO+(NUM-1 )*HDAT 
RETURN 
END 

SUBROUTINE  POINTS 
INTEGER  CODE 

COMMON/BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /MAT  P/RO ( 1 2 ) ,E(9,12) ,EE(9) ,ETA( 12) 

COMMON /ARG/XXX(  10)  ,mO0)  ,S(24)  ,XX(3)  ,1^(3) , 

1CRZ(6,6) ,XI(10) ,SIG(12) ,N,M 
COMMON /NPDATA/X( 200) ,Y(200) ,CODE(200) ,NPNUM( 10 ,20) 

COMMON /ELDATA/BETA (12) ,ALPHA(12) ,TH(12) ,IX(200,4) ,MATRIL(12) 
COMMON/RESULT/D (6, 6 ) ,C(6,6) ,CNS (6, 6) 

COMMON /TD /IMIN (100), IMAX (100), JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM,TSUM( 6 ) , SIGMA ( 12,50,6) ,D3IG(6 ) ,F(600 ) 

COMMON /DISP/DELN1 (600) ,DELN(600) , DEL (600) ,GNM1 (600 ) ,GNM2(600 ) 
COMMON/MASS/A(600) ,B(600) ,CM(8) 

COMMON /PLYLD /SIGY (12,6), DEPS (6 ) , EPS ( 1 2 , 50 , 6 ) 

DIMENSION  AX( 10 ,20) , AY( 1 0 ,20) ,BLKANG ( 1 2) ,BLKALF( 12) 
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EQUIVALENCE  (X( 1 ) , AX) , ( Y( 1 ) , AY) 

*********************************** 

ESTABLISH  NODAL  POINT  INFORMATION 

*********************************** 

NEL=0 
N0DSUM=0 
DO  100  J=1 ,MAXJ 
NSTART=IMIN(J) 

NSTOP=IMAX(J) 

DO  100  I=NSTART ,NSTOP 
100  N0DSUM=N0DSUM+1 
NELSUM=0 
JJMAX=MAXJ-1 
DO  110  JJ=1 , JJMAX 
NS'IOP=MINO(IMAX( JJ )  ,IMAX(JJ+1 ) )-1 
NSTART=MAXO (IMIN ( JJ ) , IMIN ( JJ+1 ) ) 

DO  110  II=NSTART ,NSTOP 
110  NELSUM=NELSUM+1 
NUMNP=NODSUM 
NUMEL=NELSUM 
DO  120  J=1,MAXJ 
NSTART=IMIN ( J) 

NSTOP=IMAX(J) 

DO  120  I=NSTART , NSTOP 
NPNJM(I , J)=NODE (1 , J) 

NP=NPNUM(I ,  J) 

X(NP )=AX(I , J) 

120  Y(NP )=AY(I , J ) 

c*  *********************************** 

C  READ  AND  ASSIGN  BOUNDARY  CONDITIONS 

Q*  *********************************** 

C  INITIALIZE 

C*  *********************************** 

NUMN= (NUMLA+1 )»NUMNP 
DO  130  1=1 ,NUMN 
CODE(I)=0 
130  CONTINUE 

IF(NBC.EQ.O)  GO  TO  210 
DO  200  IBCON=1 ,NBC 
READ(5,1002)N1 ,N2tICN 
DO  200  I=N1 ,N2 
CODE(I)=ICN 
200  CONTINUE 
210  MPRINT=0 

DO  230  J=1 ,MAXJ 
NSTART=IMIN ( J ) 

NSTOP=IMAX( J ) 

DO  230  I=NSTART, NSTOP 
NP=NPNUM(I , J) 

IF(MPRINT.NE.O)  GO  TO  220 
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WfiITE(6 ,2000 ) 

MPRINT=59 

220  MPRINT=MPRINT-1 

230  WfiITE(6 ,2001 )  I , J ,NP  , CODE (NP) ,X(NP ) ,  Y (NP ) 

DO  310  IMTL=T,NUMLA 

READ  (5,1000)  MTL, BETA ( IMTL ), ALPHA (IMTL) , ETA (IMTL) 

READ (5 ,2444) (SIGY (IMTL ,1) ,1=1 ,6) 

2444  FORMAT (6E 12. 6) 

IF (SIGY (IMTL , 2) ,NE .0 .0 )  GO  TO  244 
DO  242  1=2,3 

242  SIGY (IMTL , I)=SIGY (IMTL , 1 ) 

DO  24^  1=4,6 

243  SIGY (IMTL, I )=SIGY (IMTL, 1 )/SQRT( 3.0) 

244  CONTINUE 

310  MATRIL ( IMTL ) =MTL 

WRITE(6 ,2003)  (CODE (I) ,1=1 ,NUMN) 

2003  FORMAT (2513) 

C*  ****************************** 

C  ESTABLISH  ELEMENT  INFORMATION 

C*  ****************************** 

JJMAX=MAXJ-1 
N=0 
MTL  =  1 
KTL=1 

DO  440  JJ=1,JJMaX 

NSTOP=MINO (IMAX ( JJ ) ,IMAX( JJ+1 ) )-1 

NSTART=MAXO (IMIN ( JJ ) ,IMIN (JJ+1 ) ) 

DO  440  II=NSTART ,NSTOP 
NEL=NEL+1 

420  i=NPNUM(II , JJ ) 

J=I+1 

K=NPNUM(II+1 ,JJ+1 ) 

L=K-1 
M=NEL 
lX(M, 1 )=I 
IX(M,2)=J 
IX(M,3)=K 
IX(M,4)=L 
440  CONTINUE 

IF(NUMNP .GT .2000 )  WRITE(6 ,2002) 

1000  FORMAT  (  15 , 2F10 . 0 ,E12 . 6 ) 

1002  FORMAT (21 5, 110) 

2000  FORMAT  (59H1  I  J  NP  TYPE  X-ORDINATE 

1TE  ) 

2001  FORMAT  (215 , 16 , 112 ,F1 3 .6 ,F1 4 . 6 , £26 .7 ,E24 . 7 ,E24 .7 ) 

2002  FORMAT  (35H  BAD  INPUT  -  TOO  MANY  NODAL  POINTS) 

RETURN 

END 

SUBROUTINE  QUAD 
INTEGER  CODE 


*  *  *  * 

*  *  *  * 


Z-ORDINA 
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REAL  NUSN ,NUTN , NUTS , NUNS  ,NUNT ,NUST 
DIMENSION  DUMMY (6, 6) , DUMMY 1(6,6) 

COMMON/BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON / MAT  P / RO ( 1 2 ) ,E(9,12) ,EE(9) ,ETA( 12) 

COMMON / ARG/XXX (10), YYY (10) ,S(24),XX(3) ,YY(3) , 

1CRZ(6 ,6) ,XI(10) ,SIG(12)  ,N,M 
COMMON /N PD AT A /X (200 ) , Y (200 ) , CODE ( 200 ) ,NPNUM( 10,20) 

COMMON /ELDATA/BETA( 12) ,ALPHA(12) ,TH(12) ,IX(200 ,4) ,MATRIL(12) 

COMMON /RESULT /D( 6, 6) ,C(6 ,6) ,CNS(6 ,6) 

COMMON /TD/IMINC 100 ) ,IMAX( 100), JMIN (25) ,JMAX(25) , MAXI ,MAXJ ,NMTL , NBC 
COMMON /BASIC2 /BET ,H,HN,HH(3),HT(2), TIME , NLA Y , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM , TSUM ( 6 ) , SIGMA ( 1 2 , 50 , 6 ) , DSIG ( 6 )  , F ( 60 0 ) 

COMMON /DISP/DELN 1 (600 ) ,DELN (600) ,DEL(600) ,GNM1(600) ,GNM2(600) 
COMMON /MASS/A( 600) ,B(600) ,CM(8) 

IF(NTIME.EQ.I)  GO  TO  10 
READ(1)  ((CRZ(I,J) ,J=1 ,6) ,1=1 ,6) 

GO  TO  151 
10  CONTINUE 
II =IX(N , 1 ) 

J1=IX(N ,2) 

K1 =IX(N ,3) 

L1=IX(N  ,4) 

MTYPE=MATRIL(NLAY) 

C*  *********************************** 

C  INTERPOLATE  MATERIAL  PROPERTIES 

q*  *********************************** 

DO  100  1=1,9 
100  EE(I)=E(I ,MTYPE) 

DO  110  1=1,6 
DO  110  J=1 ,6 
CNS(I , J)=0.00 
C(I,J)=0.00 
110  D(I , J )=0 .00 

C*  *********************************** 

C  FORM  STRESS-STRAIN  RELATIONSHIP  IN  N-S-T  SYSTEM 

q*  *********************************** 

NUNS=EE(4) 

NUNT=EE(5) 

NUST=EE (6 ) 

NUSN=(EE(2)*NUNS)/EE(1 ) 

NUTN=(EE(3)*NUNT )/EE( 1 ) 

NUTS=(EE(i)*NUST)/ES(2) 

DIV=1 . 00-NUNS*NUSN-NU3T*NUTS-NUNT*NUTN-NUSN*NUNT*NUTS 
1-NUNS*NUTN*NUST 

CNS ( 1 , 1 ) =EE ( 1 ) * ( 1 . 00-NUST*NUTS ) /DI V 
CNS ( 1 , 2) =EE(2 )* (NUNS+NUNT*NUT3) /DIV 
CNS ( 1 , 3 ) =EE ( 3 ) * ( NUNT+NUNS*NUST ) /DIV 
CNS(2 , 1 )=CNS( 1 ,2) 

CNS ( 2 , 2 ) =£E ( 2 ) * ( 1 . 00 -NUNT*NUTN ) /DIV 
CNS ( 2 , 3 ) =EE ( 3 ) * (NUST+NUSN  *NUNT ) /DIV 
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CNS(3,1)=CNS(1,i) 

CNS(3,2)=CNS(2 ,3) 

CNS ( 3 , 3 ) =EE ( 3 ) * ( 1 . 00-NUNS*NUSN ) /DIV 
CNS (4 , 4 )=EE(7 ) 

CNS (5 ,5)=EE(8) 

CNS (6 ,6)=EE(9) 

C  SET  UP  STRAIN  TRANSFORM  TO  N-S-T  SYSTEM 
3INA=SIN (ALPHA(M) ) 

COSA=COS(ALPHA(M) ) 

S2=SINA**2 
C2=COSA**2 
SC=SINA*C05A 
D( 1 , 1 )=C2 
D(1,3)=S2 
D( 1 ,6)=-SC 
D(2 , 1 )=S2 
D(2,3)=C2 
D(2 ,6)=SC 
D(3,2)=1 .00 
D(4 , 1 )=2.00*SC 
D (4 , 3 )=-2 . 00*3C 
D(4 ,6)=C2-S2 
D (5 ,4)=SINA 
D(5,5)=COSA 
D(6,4)=COSA 
D(6,5)=-SINA 

C  SET  UP  STRAIN  TRANSFORMATION  TO  R-Z-T  SYSTEM 
3INB=SIN (BETA(M) ) 

C03B=C0S(BETA(M) ) 

S2=SINB»»2 

C2=C0SB**2 

SC=SINB*COSB 

C(1 ,1)=S2 

C(1 ,2)=C 2 

C(1 , 4 )=SC 

C(2 , 1 )=C2 

C(2,2)=S2 

C(2,4)=-SC 

C(3,3)=1.00 

C(4,1)=-2.00*SC 

C(4,2)=2.00*SC 

C(4 ,4)=S2-C2 

C(5,5)=SINB 

C(5,6)=-COSB 

C(6,5)=COSB 

C(6,6)=SINB 

C  CALCULATE  CRZ  MATRIX 
DO  120  1  =  1  ,6 
DO  120  J  =  1  ,6 
DUMMY(I , J)=0 .00 
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DO  120  K=1 ,6 

120  DUMMY (I , J)=DUMMY(I , J)+D(I ,K)*C(K , J) 

DO  130  1=1,6 
DO  130  J=1 ,6 
DUMMY 1 (I,J)=0.00 
DO  130  K=1 ,6 

130  DUMMY 1 (I, J)=DUMMY' (I , J)+CNS(I ,K)*DUMMY(K , J) 

DO  140  1=1,6 
DO  140  J=1 ,6 
DUMMY ( I , J ) =0 . 00 
DO  140  K=1  ,6 

140  DUMMY ( I , J ) =DUMMY ( I , J)+D(K ,I)*DUMMY1 (K,J) 

DO  150  1=1 ,6 
DO  150  J=1 ,6 
CRZ(I,J)=0.00 
DO  150  K=1 ,6 

150  CRZ(I , J)=CRZ(I , J)+C(K ,I)*DUMMY (K , J) 

WRITE(1 )  ( (CRZ(I , J) ,J=1 ,6), 1=1 ,6) 

151  CONTINUE 

DO  200  1=1,4 
MM=  IX(N,I) 

XXX(I)=  X(MM) 

XXX (1+4) =  X(MM) 

YYY (I )=  Y(MM) 

200  YYY(I+4)=  Y (MM) 

XXX ( 9 ) =  (XXX(1)+  XXX(2 )+  XXX(3)+  XXX(4))/4.00 
YYY(9)=  (YYY ( 1 )+  YYY(2)+  YYY(3)+  YYY(4))/4.00 
XXX(10)=  XXX(9 ) 

YYY (10)=  YYY(9) 

DO  250  1=1,8 
CM(I)=0 .00 
250  CONTINUE 

DO  900  IA=1 ,24 
900  S(IA)=  0.00 

CALL  TRISTF( 1,2) 

CALL  TRISTF(2,3) 

CALL  TRISTF(3 ,4) 

CALL  TRISTF(4 , 1 ) 

RETURN 

END 

SUBROUTINE  STIFF 
INTEGER  CODE 

COMMON / BASIC /VOL , NL  MNP , NUMEL , NUMLA , NCG 
C0MM0N/MATP/R0(12) ,E(9,12) ,EE(9) ,ETA( 12) 
COMMON/ARG/XXX(10),YYY(10),S(24),XX(3),YY(3), 

1CRZ(6 ,6) , XI ( 1 0 ) , SIG ( 1 2 ) , N , M 

COMMON/NPDATA/X ( 200 ) , Y ( 200 ) , CODE ( 200 ) , NPNUM( 1 0 , 20 ) 
COMMON/ELDATA/BETA( 12 ) , ALPHA ( 12), TH (12), IX (200 ,4) ,MATRIL( 1 2) 
COMMON/RESULT/D (6 ,6) ,C(6 ,6) ,CNS(6,6) 
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COMMON /TD/1MIN (100) ,IMAX(100) ,JMIN (25) ,JMAX(25) , MAXI, MAXJ,NMTL, NBC 
COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLA I , IFLAG , IT , NTIME 
COMMON /SIGM/BSUM,TSUM( 6) ,SIGMA( 12,50 ,6) ,DSlG(6) ,F(600) 

COMMON /DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 

COMMON /MASS/ A (600) ,B(600) ,CM(8) 

COMMON/MASS 1/XMINV( 600) ,EEKM(1 3) ,F1 (600) ,F2(600) 

IFLAG=1 
DO  50  1  =  1  ,  IT 
GNM2(I ) =GNM1 (I ) 

GNM1 (I )=B( I ) 

50  CONTINUE 

DO  100  1=1 , IT 
B(I)=  0.00 
100  CONTINUE 
REWIND  1 
REWIND  2 
REWIND  3 

DO  340  N=1 , NUMEL 
DO  340  M=1,NUMLA 
NLAY=M 
CALL  QUAD 
DO  340  1  =  1  ,4 

II  =  3*IX(N,I)+  3*(M-1)*NUMNP 

J=  I 1-2 

DO  340  K=J , II 

JJ=K-II+3*I 

B(K)=  B(tC)+S(JJ) 

KK=  K+3*NUMNP 
B(KK)=  B(KK)+S(JJ+12 ) 

IF ( NTiME. GT . 1 )GO  TO  340 
A(K)=A(K)+CM(I) 

A(KK)=A(KX)+CM(l+4) 

340  CONTINUE 

DO  4001=1, IT 
400  B(I)=-B(I) 

RETURN 

END 

SUBROUTINE  STRESS 
INTEGER  CODE 

COMMON /BASIC/ VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/NPR/NPRINT 

COMMON /MATP/RO( 12) ,E(9, 12) ,E£(9) ,ETA(12) 

COMMON /ARG/XXX( 10) ,YYY ( 10 )  ,S(24) ,XX(3) , YY ( 3 ) , 

1CRZ(6 ,6) ,XI( 10 ) ,SlG( 12 )  ,N ,M 

COMMON/NPDATA/X(200) ,Y(200) , CODE (200) ,NPNUM( 10 , 20 ) 
COMMON/ELDATA/BETA( 12) , ALPHA ( 12) ,TH( 1 2) , IX (200 , 4 ) ,MATR1L (12) 
COMMON/RESULT/D (6 ,6) ,C(6 ,6) ,CNS(6 , 6) 

COMMON /TD/IMIN ( 100 ) , IMAX( 100 ) , JMIN (25) , JMAX(25 ) ,MAXI ,MAXJ ,NMTL ,NBC 
COMMON/BASIC2/BET , H ,HN , HH( 3 ) ,HT (2 ) , TIME ,NLAY , IFLAG , IT , NTIME 
C0MM0N/3IGM/BSUM,T3UM(6) ,SIGMA( 12,50,6) ,DSIG(6) ,F(600) 
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C0MM0N/DISP/DELN1 (600) ,DELN(600) , DEL (600) ,GNM1 (600) ,GNM2(600) 
C0MM0N/MASS/A(600) ,B(600) ,CM(8) 

COMMON /DELT/XDEL (4 ,  18 ) 

COMMON/DELTRI/DELTA ( 4 , 1 8 ) 

COMMON/PLYLD/SIGY (12,6), DEPS ( 6 ) , EPS ( 1 2 , 50 , 6 ) 

COMMON/SIGZ/DSIGZ ( 1 2 ) 

DIMENSION  DDEL(600) 

C  CALCULATE  CHANGES  IN  DISPLACEMENTS 
DO  100  1=1, IT 
ddel(i)=del(i)-deln (I) 

100  CONTINUE 

DO  101  1=1,4 
DO  101  J  =  1  , 1 8 
XDEL(I , J)=0.00 
DELTA(I , J)=0.00 

101  CONTINUE 
IF LAG =2 
REWIND  1 
REWIND  2 
REWIND  3 

C  WRITE(6 ,102) 

102  FORMATC'I  THESE  ARE  THE  STRESSES  CALCULATED  FROM  THE  STRAINS  AND 

1  PLASTIC  YIELD  WITHOUT  MODIFICATION"/ 

2  "  LAYER ",2X," ELEMENT ",5X," SIGMA  X  ",6X," SIGMA  Y  ",6X, 

3  "SIGMA  Z  ",6X, "SIGMA  XY ", 6X ," SIGMA  YZ" , 6X , "SIGMA  XZ") 

DO  200  N=1 ,NUMEL 

CALL  ZSTRESS 
DO  200  M=1 ,NUMLA 
NLAY=M 
DO  150  1=1,4 
11=1+1 

IF(I.EQ.4)  11=1 
MM=IX(N ,1) 

MM1 =IX(N ,11 ) 

II=3*MM+3*(M-1 )*NJMNP 

II1=3*MM1+3*(M-1 )*NUMNP 

J=II-2 

J 1=11 1-2 

DO  140  K=J , II 

IK=K-J+1 

XDEL (I , IE ) =DDEL (K ) 

DELTA ( I , IK)=DELN (K) 

KK=K+3*NUMNP 

XDEL ( I , IK+6 ) =DDEL ( KK ) 

DELTA ( I , IK+6 ) =DELN ( KK ) 

140  CONTINUE 

DO  145  K=J1 ,111 
IK=K-J1+1 

XDEL ( I , IK+3 ) =DDEL ( K ) 

DELTA ( I , IK+3 ) =DELN ( K ) 
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KK=K+3*NUMNP 

XDEL(I,IK+9)=DDEl(KK) 

DELTA ( I , IK+9 ) =DELN ( KK ) 

145  CONTINUE 
150  CONTINUE 

C  FIND  DISPLACEMENT,  BY  AVERAGING,  AT  CENTER  OF  EACH  ELEMENT 

XDEL(1 ,13)=(XDEL(1 ,1 )+XDEL(2,1)+XDEL(3,1 )+XDEL(4,1))/4.00 
XDEL(1 , 14 )= (XDELC 1 ,2)+XDEL(2,2)+XDEL(3,2)+XDEL(4,2))/4.00 
XDEL(1 ,15)=(XDEL(1 ,3)+XDEL(2,3)+XDEL(3,3)+XDEL(4,3))/4.00 
XDEL(1 ,16)=(XDEL(1 ,7)+XDEL(2,7)+XDEL(3,7)+XDEL(4,7))/4.00 
XDEL(1 ,17)=(XDEL(1 ,8)+XD£L(2,8)+XDEL(3,8)+XDEL(4,8))/4.00 
XDEL (1 , 1 8 ) = (XDEL ( 1 , 9 )+XDEL(2 , 9 )+XDEL(3 , 9 )+XDEL(4 , 9 ) ) /4 . 00 
DELTA (1 , 1 3 )= (DELTA (1 , 1 ) +DELTA (2,1) +DELTA (3,1) +DELTA (4,1)) 
1/4.00 

DELTA(1 , 14)  =  (DELTA(1 ,2)+DELTA(2,2)+DELTA(3 ,2)+DELTA(4,2) ) 
1/4.00 

DELTA ( 1 , 1 5 ) = ( DELTA (1,3) +DELT A ( 2,3) +DELTA (3,3 )+DELT A ( 4,3)) 
1/4.00 

DELTA(1 , 16)=(DELTA(1 ,7 )+DELTA(2 , 7 )+DELTA(3 ,7 )+DELTA(4 ,7 ) ) 
1/4.00 

DELTA(1 , 17)=(DELTA(1 , 8 )+DELTA(2 , 8 )+DELTA( 3 , 8 )+DELTA(4 , 8 ) ) 
1/4.00 

DELTA ( 1 , 1 8 ) = ( DELTA ( 1 , 9 ) +DELTA (2,9 )+DELTA (3,9) +DELT A ( 4,9)) 
1/4.00 

C  CALCULATE  AVERAGE  STRAIN 
DO  175  1=2,4 
DO  175  J= 1 3 , 1 8 
XDEL ( I , J ) =XD£L ( 1 , J ) 

DELTA (I , J)=DELTA( 1 , J) 

175  CONTINUE 
179  BSUM=0 . 00 

DO  180  1=1 ,6 
TSUM(I)=0.00 
180  CONTINUE 
CALL  QUAD 
DO  190  1=1 ,6 
D£P3(I)=TSUM(I)/BSUM 
190  EPS(M,N , I)=EPS(M,N , I )+DEPS(I ) 

999  CONTINUE 

C  CALCULATE  STRESS  INCREMENT 
DO  195  1=1 ,6 
DSIG(I)=0.00 
DO  195  J=1 ,6 

195  DSIG(I)=  DSIG (I )+CRZ(I , J)*D£PS( J ) 

CALL  YIELD 
200  CONTINUE 
RETURN 
END 

SUBROUTINE  TRISTF(II , JJ) 

INTEGER  CODE 
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COMMON /BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON/MATP/RO(12),E(9,12) ,EE(9) ,ETA( 12) 

COMMON/ARG/XXX(10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 ,6)  f XX ( 1.0 )  ,SIG(12) ,N,M 
COMMON /NPDATA/X’(200)  ,Y(200 ) , CODE (200)  ,NPNUM(  10 ,20) 

COMMON /ELDATA/BET A ( 1 2 ) , ALPHA ( 1 2 ) , TH ( 1 2 ) , IX ( 200 , 4 ) , MAT  RIL ( 1 2 ) 

COMMON /TD /IMIN (100), IMAX ( 1 00 ) , JMIN ( 25 ) , JMAX ( 25 ) , MAXI , MAX J , NMTL , NBC 
COMMON /BASI C2 /BET , H , HN , HH ( 3 ) , HT ( 2 ) , T IME , NLA Y , IFLAG , IT , NT IME 
COMMON /SIGM/BSUM , TSUM ( 6 ) , SIGMA (12,50,6), DSIG (6) ,F(600) 

COMMON /DISP/DELN 1(600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
COMMON /MASS/ A (600) ,B(600) ,CM(8) 

COMMON /DELT /XDEL (4,18) 

COMMON /DELTRI /DELTA (4,18) 

DIMENSION  Q(6 , 1 8 ) , C ( 1 8 , 1 8 ) ,XY(9,9) , SAVE (18, 3) ,BT(18,6) ,BTSIG(18) , 
1CBAR( 10 , 18)  ,Z( 10 , 10) ,  CBD(10),CC(10,4) ,BLT(18,4) , 

2BADD (18,6) 

DIMENSION  BLTSIG(18) ,BSIG(18) 

DIMENSION  DUM( 10,10) 

NLA=M 

NEL=N 

DO  50  1=1 ,6 

50  SIG ( I ) =SIGMA ( M , N , I ) 

THICK=TH(NLAY) 

ISUB=  9 

XX(1)=  XXX(II) 

XX(2)=  XXX(JJ) 

XX ( 3 ) =  XXX(ISUB) 

YY ( 1 ) =  YYY(II) 

YY(2 )=  YYY(JJ) 

YY  ( 3 ) =  YYY(ISUB) 

IF(NTIME.EQ.I)  GO  TO  40 
READ (2)  ((BT(I,J) ,J=1 ,6) ,1=1 ,18) 

GO  TO  801 
40  CONTINUE 
100  CALL  INTER 
150  CONTINUE 

DO  300  IA=1 ,6 
DO  300  JA=1 , 18 
300  Q(IA, JA)=0 .00 

Q(1,2)=  XI ( 1 ) *  THICK 

Q ( 1  ,11)=  XI(1  )*( ( THICK )**2)/2. 00 

Q(3 , 18 )=  XI(3)»  THICK 

Q(5 , 14)=  XI(2)»  THICK 

Q(2 ,6)=  Q(1 ,2) 

Q(3 , 16 )=  Q( 1 ,2) 

Q(3 , 17 )=  Q(5 , 14) 

Q(4,3)=  Q( 1 ,2) 

Q(4 ,5)=  Q(1 ,2) 

Q(5,9)=  Q(1 ,2) 

Q(5,13)=  Q(1,2) 
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Q(6 ,8)=  Q( 1 , 2) 

Q(6 , 10 ) =  Q(1 ,2) 

Q(2 , 15)=  Q(1 , 11 ) 

Q(4,12)=  Q ( 1 ,11) 

Q(4 , 14 )=  Q ( 1 ,11) 

Q(5,18)=  Q ( 1 ,11) 

Q(6 , 17 )=  Q( 1 , 1 1 ) 

Q(5 , 15 )=  Q(5 , 18 ) 

Q(6 , 12)=  Q  C  3 » 1 8 ) 

Q(6 , 1 1 )=  Q ( 5 , 1 4 ) 

THE  FOLLOWING  EXPRESSION  IS  ACTUALLY  TWICE  THE  AREA  OF  THE  TRIANGLE. 

AREA=  XX( 1 )*YY (2)-XX(2 )*YY ( 1 )+XX(2)*YY(3 )-XX(3 )*YY (2)+XX(3 )*YY ( 1 )- 
1XX( 1 )*YY(3) 

ZERO  C  MATRIX 

DO  400  1=1,18 
DO  400  J  =  1 ,18 
400  C(I,J)=  0.00 
FAC=1 /AREA 
M=  1 
N=  1 
1=  II 
J=  JJ 
GO  TO  250 
75  CONTINUE 

FAC= 1 / ( AREA*THICK ) 

M=  10 
N  =  10 
1=11+4 
J=JJ+4 


FILL  C  INVERSE  BY  PARTS 

250  C(M,N)=  (XXX(J)*  YYY(ISUB)-  XXX(ISUB)*  YYY(J))*FAC 
C(M,N+3)=  (XXX(ISUB)*  YYY (I )-  XXX(I)*  YYY(ISUB) )*FAC 
C(M,N+6)=  (XXX(I)*  YYY(J)-  XXX(J)*  YYY(I))*FAC 
C(M+1 ,N)=(YYY(J)-  YYY(ISUB))*FAC 
C(M+1 , N+3 )= (YYY (iSUB)-  YYY(I))*FAC 
C(M+1,N+6)=(YYY(I)-  YYY(J))*FAC 
C(M+2 ,N)= (XXX (ISUB)-  XXX(J))*FAC 
C(M+2 ,N+3 )= (XXX(I)-  XXX ( ISUB ))*FAC 
C(M+2 ,N+6 )= (XXX( J )-  XXX (I) )*FAC 
C(M+3 ,N+1 )=  C(M,N) 

C(M+3,N+4)=  C(M,N+3) 

C(M+3 ,N+7 )=  C(M,N+6 ) 

C(M+4,N+1 )=  C(M+1 ,N) 

C(M+4,N+4)=  C(M+1 ,N+3 ) 
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C(M+4 ,N+7 )=  C(M+1 ,N+6) 
C(M+5,N+1)=  C ( M+2 , N ) 
C(M+5 ,N+4)=  C(M+2,N+3) 
C(M+5 ,N+7 )=  C(^+2,N+6) 
C(M+6 ,N+2)=  C(M,N) 

C ( M+6 , N+5 ) =  C(M,N+3) 

C ( M+6 , N+8 ) =  C(M,N+6 ) 

C ( M+7 , N+2 ) =  C ( M+ 1 , N ) 

C ( M+7 , N+5 ) =  C(M+1 ,N+3) 

C ( M+7 , N+8 ) =  C(M+1 ,N+6 ) 
C(M+8 ,N+2)=  C ( M+2 , N ) 

C ( M+8 , N+5 ) =  C(M+2,N+3) 
C(M+8,N+8)=  C(M+2,N+6) 
IF(ISUB.EQ.IO)  GO  TO  200 
ISUB=  10 
GO  TO  75 
200  M=NLA 
N=NEL 


C  COMPUTING  THE  LOWER  LEFT  QUADRANT  OF  THE  C  MATRIX 
C 

DO  600  1=1,9 
DO  600  J=1 ,9 
IJ=I+9 

600  C(IJ , J)=-C(I , J)/THICK 
C 

C  REARRANGE  THE  DISPLACEMENT  MATRIX 
C 

DO  730  IA=1 , 18 
DO  730  JA=1 ,3 
JB=  JA+  6 

730  SAVE(IA , JA)=  C(IA,JB) 

DO  740  IA=1 ,18 
DO  740  JA=7 , 12 
JB=JA+3 

740  C(IA, JA)=  C(IA, JB) 

DO  750  IA=1 , 18 
DO  750  JA=1 ,3 
JB=JA+12 

750  C(1A,JB)=  SAVE(IA, JA) 

C 

C  CALCULATE  THE  BO  TRANSPOSED  MATRIX 
C 

DO  800  JA=1 , 18 
DO  800  IA=1 ,6 
BT(JA,IA)=0.00 
DO  800  MA=1 , 18 

800  BT ( JA, IA)=BT ( JA, IA)+C(MA , JA)*Q(IA,MA) 

WRITE(2)((BT(I,J),J=1 ,6), 1=1, 18) 

80 1  CONTINUE 
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CALCULATE  THE  BOT  *  STRESS 

DO  810  IA=1 ,18 
BTSIG(IA)=0 .00 " 

DO  810  NA=1 ,6 

8 1 0  BTSIG ( I A ) =BTSIG ( IA ) +BT ( 1A , NA ) *SIG (NA ) 

C  ADD  THE  NONLINEAR  TERMS 

C 

C 

C  CALCULATE  THE  CBAR  MATRIX 
C 

IF(NTIME.EQ.I)  GO  TO  815 

C  LN  IS  THE  NUMBER  IN  TIME  STEPS  BETWEEN  CALCULATION  FOR  LARGE  DEFORMATION 
LN=2 

IF( (NTIME/LN  )*LN .EQ.  NTIME)  GO  TO  815 
READ(3)  ( ( BLT ( I , J ) , J  =  1 ,4) ,1=1 ,18) 

GO  TO  861 
815  CONTINUE 

DO  820  IA=1 , 10 
DO  820  JA=1 , 18 
IB=  IA+7 

IF ( IA . EQ . 9 ) IB= 1 7 
IF(IA.EQ. 10)IB=18 
820  CBAR(IA,JA)=  C(IB,JA) 

C  CALCULATE  THE  ZBAR  MATRIX 
C 

DO  830  IA=1 , 10 
DO  830  JAsI ,10 
830  Z(IA, JA)s  0.00 

Z(1  ,1)=  XI ( 1 )*THICK 

Z(2,2)=  Z(1,1) 

Z(3,3)=  Z(1,1) 

Z(6 ,6)=  Z(1  ,1) 

Z(3,4)=  XI (2 )*THICK 
Z(4 , 3)=  Z(3 ,4) 

Z(6 ,7 )=  Z(3,4) 

Z(7 ,6)=  Z(3 ,4) 

Z(9 ,9)=  (THICK**3)/3*XI(1 ) 

Z(10,10)=  Z(9 ,9) 

Z ( 1  ,9)=  (TH1CK**2)/2*XI(1 ) 

Z(2 ,10)=  Z( 1 ,9) 

Z(9 , 1 )=  Z(1 ,9) 

Z( 10 ,2)=  Z(1 ,9) 

Z(3,5)=  XI ( 3 ) *  THICK 
Z ( 5 ,3)=  Z(3 ,5) 

Z(6 ,8)=  Z(3,5) 

Z(8,6)=  Z(3,5) 
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Z(4 ,4)=  XI (5  )*  THICK 
Z ( 7 ,7)  =  Z(4,4) 

Z(5 , 5) =  XI(6)«  THICK 
Z ( 8 , 8 )  =  Z(5  >5) 

Z(4,5)=  XI (4 )*  THICK 
Z(5 ,4)=  Z(4 ,5) 

Z(7,8)=  Z(4 ,5) 

Z  ( 8 , 7 ) =  Z(4 ,5) 

C 

CALCULATING  THE  CC  MATRIX 

DO  840  IA=1 , 10 
CBD(IA)=  0.00 
DO  840  JA=1 ,18 

840  CBD(IA)=  CBD(IA)+  CBAR(IA,JA)«  DELTA(II , JA) 
DO  850  IA=1  , 10 
DO  850  JA=1  ,4 
850  CC(IA, JA)=0 .00 
CC(1 ,1)=  CBD(1) 

CC(2 ,4)=  CBD( 1 ) 

CC(1 ,4)=  CBD(2 ) 

CC(2 , 2)=  CBD(2 ) 

CC(i ,3)=  CBD ( 3 ) 

CC(4 ,3)=  CBD(4 ) 

CC(5,3)=  CBD ( 5 ) 

CC (6,3)=  CBD (6) 

CC(7,3)=  CBD(7 ) 

CC(8,3)=  CBD(8 ) 

CC (9,1)=  CBD(9 ) 

CC(10,4)=  CBD(9 ) 

CC (9,4)=  CBD(IO) 

CC( 10,2)=  CBD(IO) 

DO  855  1=1,10 
DO  855  J=1 ,4 
DUM(I , J)=0.00 
DO  855  K=1 ,10 

DUM(I,J)=DUM(I,J)+Z(I,K)*CC(K,J) 

855  CONTINUE 

DO  860  1  =  1 ,18 
DO  860  J  =  1 ,4 
BLT (I , J )=0 .00 
DO  860  K=1 ,10 

BLT (I , J)=BLT (I , J )+CBAR(K ,I)*DUM(K , J ) 

860  CONTINUE 

WRITE(3 )  ( ( BLT ( I , J ) ,J=1 ,4) ,1=1 ,18) 

861  CONTINUE 

DO  870  IA=1 ,18 
BLTSIG(IA)=  0.00 
DO  870  JA=1 ,4 

870  BLTSIG(IA)=  BLTSIG(IA)+  BLT(IA,JA)»  SIG(JA) 
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DO  880  IA= 1 , 18 

BSIG(IA)=  BTSIG(IA)+  BLTSIG(IA) 

880  CONTINUE 

C  ADDING  THE  TRIANGULAR  B3IG  MATRICES  TO  FORM  THE  QUADRILATERAL  MATRIX 
IF(IFLAG.NE.2)  GO  TO  1600 
C  EVALUATE  STRAIN  INCREMENT 

BSUM=BSUM+AREA/2 . 0*THICK 
DO  525  1=1,18 
DO  525  J=1 ,4 

BADD(I , J)=  BT(I,J)+  BLT (I , J ) 

525  CONTINUE 

DO  526  1=1 , 18 
DO  526  J=5 ,6 

526  BADD ( I , J ) =  BT(I , J) 

DO  550  1=1 ,6 

DO  550  J  =  1 , 18 

550  TSUM(I ) =T3UM( I )+BADD ( J , I ) *XDEL( II , J ) 

GO  TO  1000 
1600  CONTINUE 
K=  3*H-2 
L=  3*JJ-2 
DO  910  MA=1 ,2 
MJ=  (MA-1 )*  6 
NJ=  (MA-1)*  12 
DO  920  IA=1  ,3 
JA=  IA-1 

S(K+  NJ+  JA) =  S(K+  NJ+  JA)+  BSIG(MJ+  IA) 

920  S(L+  NJ+  JA)=  S(L+  NJ+  JA)+  BSIG(MJ+  IA+  3) 

DO  910  IA=1 , 3 
DO  910  JA=1,4 
IJ=  (JA-1)*3 

910  S(NJ+IJ+IA)=  S(NJ+IJ+IA)+  BSIG(12+(MA-1)*3+IA)/4.0 
IF(NTIME.GT.I)  GO  TO  1000 
C  ASSEMBLE  ELEMENT  MASS  MATRIX 

DEN=RO(MATRIL(NLAY))*THICK 

XMASS=AREA*DEN/8.00 

CM ( 1 1 ) =  CM (I I ) +XMASS 

CM( JJ )=CM( JJ )+XMASS 

CM(II+4)=CM(II) 

CM( JJ+4 ) =CM( JJ ) 

1000  CONTINUE 
RETURN 
END 

SUBROUTINE  YIELD 
COMPLEX  LAMBDA,  DISC,  DISQ 

COMMON /ARG/XXX( 10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1CRZ(6 , 6) ,XI ( 10) ,SIG( 12) ,N ,M 
COMMON /SIGM/BSUM,TSUM( 6 )  ,3IGMA(  12,50,6)  ,DSIG(6)  ,F(600) 
C0MM0N/BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IFLAu , IT , NTIME 
COMMON /PLYLD/SIGY (12,6) ,DEPS(6 ) ,EPS( 12,50,6) 
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COMMON /MAT P/RO( 12 ) ,E(9 , 12) ,EE(9) ,ETA( 12) 

COMMON/SIGZ/DSIGZC 12) 

DIMENSION  SIGT(6),SIGYB(3) ,TSIG(6) ,TSIGB(6) 
SIG3=SIGMA(M,N,3)+DSIGZ(M) 

C  FIND  COMBINATION  YIELD  STRESSES  (SIGYB) 

SIGYB(1 )=1 .0/SIGY (M, 1 )**2-1 .0/SIGY(M,2)»»2-1 .0/SIGY(M,3)**2 
SIGYB(2)=1 .O/SIGY (M,2)**2-1 .O/SIGY (M, 1 )«»2-1 .0/SIGY(M,3)**2 
SIGYB(3)=1 .O/SIGY (M , 3 ) **2—1 .0/SIGY(M, 1 )»*2-1 .0/SIGY(M,2)»»2 
C  GET  TEST  STRESSES  FROM  ELASTIC  ANALYSIS  (SIGT) 

DO  15  1  =  1  ,6 

15  SIGT(I)=SIGMA(M,N,I)+DSIG(I) 

C  TEST  FOR  YIELDING  BY  HILL'S  CRITERION 
YLD=0 .0 
DO  20  1=1 ,6 

20  YLD=YLD+SIGT ( I ) *»2/SIGY (M , I ) **2 

YLD= YLD+SIGYB ( 1 ) *SIGT ( 2 ) »SIGT ( 3 ) +SIGYB ( 2 ) «3IGT ( 1 ) «SIGT ( 3 ) 

1  +SIGYB ( 3 ) *SIGT ( 1 ) *SIGT ( 2 ) 

IF(YLD .LT  .  1 .0 .OR .YLD.EQ. 1.0)  GO  TO  390 
C  CORRECT  STRESSES  FOR  PLASTICITY 

C  FIND  TEST  STRAIN  INCREMENT  VALUES  (TSIG  *  LAMBDA  =  DEPS  PLASTIC) 
TSIG(1 )=SIGT(1 )/SIGY(M,1 )**2+(SIGYB(3)*SIGT(2)+SIGYB(2)« 

1  SIGT (3 ) ) /2 .0 

TSIG ( 2 ) =SIGT ( 2 ) /SIGY ( M , 2 ) * *2+ ( SIGYB ( 3 ) *SIGT (1 ) +SIGYB (1 ) » 

1  SIGT (3 ) ) /2 .0 

TSIG ( 3 ) =SIGT ( 3 ) /SIGY ( M , 3 ) **2+ ( SIGYB ( 2 ) »SIGT ( 1 ) +SIGYB (1 ) * 

1  SIGT (2 ) ) /2 . 0 
DO  25  1=4,6 

25  TSIG ( I ) =SIGT ( I ) / ( SIGY ( M , I ) **2 ) 

DO  35  1=1,6 
TSIGB(I)=0.0 
DO  35  J=1 ,6 

35  TSIGB ( I ) =TSIGB ( I )+CRZ (I , J) «TSIG ( J ) 

C  FIND  A,B,C  OF  THE  QUADRADIC  EQUATION  TO  CALCULATE  LAMBDA 
AQ=0.0 
DO  45  1=1,6 

45  AQ=AQ+TSIGB( I )**2/SIGY (M , I ) ««2 
AQ=AQ+SIGYB ( 1 ) »TSIGB ( 2 ) «TSIGB ( 3 ) 

1  +SIGYB(2)»TSIGB(1 )*TSIGB(3) 

2  +SIG YB ( 3 ) *TSIGB ( 1 ) * 

3  TSIGB(2) 

BQ=0 .0 

DO  55  1=1 ,6 

55  BQ=BQ+TSIGB(I )*3IGT (I ) /SIGY (M,I)**2 

BQ=BQ+ ( SIGYB (1 ) » ( SIGT ( 2 ) «TSIGB ( 3 ) +SIGT ( 3 ) 

1  *TSIGB ( 2 ) ) +SIGYB ( 2 ) « ( SIGT (1 ) »TSIGB ( 3 ) 

2  +SIGT ( 3 ) *TSIGB (1 ) ) + 

3  SIG YB ( 3 ) * ( SIGT (1 ) «TSIGB ( 2 ) +SIGT ( 2 ) «TSIGB ( 1 ) 

4  ))/2.0 
CQ=YLD-1 .0 
DISC=BQ»»2-AQ»CQ 
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DISQ=CSQRT (DISC) 

LAMBDA=CQ/ ( BQ+DISQ ) 

C  CALCULATE  ELASTIC-PLASTIC  STRESSES 
DO  65  1=1  ,6 

SIGMA ( M ,  N ,  I ) =SIGMA ( M ,  N ,  I ) +DSIG (I)-TSIGB(I) *LAMBDA 
65  CONTINUE 
GO  TO  400 
390  DO  395  1=1,6 
395  SIGMA ( M , N , I ) =SIGT ( I ) 

GO  TO  600 
400  CONTINUE 

C  WRITE(6,500)M,N,(SIGMA(M,N,I) ,1=1 ,6) 

500  FORMAT (218 , 2X , 6E 1 4 . 6 ) 

SIG3=(SIG3+SIGMA(M,N,3))/2. 

SIGMA ( M , N , 2 ) =SIGMA ( M , N , 2 ) -SIGMA ( M , N , 3 ) +SIG  3 
SIGMA (M,N,1 )=SIGMA(M,N,1 )-SIGMA(M,N,3)+SIG3 
3IGMA(M,N,3)=SIG3 
600  CONTINUE 
RETURN 
END 

SUBROUTINE  ZSTRESS 
INTEGER  CODE 

COMMON/BASIC/VOL , NUMNP , NUMEL , NUMLA , NCG 
COMMON /MAT P/RO (12) ,E(9,12) ,EE(9) ,ETA(12) 

COMMON /ARG/XXX( 10) ,YYY ( 10 ) ,S(24 ) ,XX(3 ) ,YY (3 ) , 

1 CRZ (6, 6), XI (10), SIG ( 1 2 ) , N , M 
COMMON /NPDATA/X( 200) ,Y(200) ,C0DE(200 ) ,NPNUM( 10 ,20 ) 
COMMON/ELDATA/BETA( 12) ,ALPHA(12) ,TH(12) ,IX(200 ,4) ,MATRIL(12) 
COMMON/TD/IMIN (100) ,IMAX(100) ,JMIN(25) , JMAX(25 ) , MAXI ,MAXJ ,NMTL , NBC 
C0MM0N/BASIC2/BET ,H,HN,HH(3),HT(2), TIME , NLAY , IFLAG , IT , NT IMS 
C0MM0N/3IGM/BSUM,TSUM(6) ,SIGMA(12,50,6) ,DSIG(5) ,F(600) 

COMMON /DISP/DELN1 (600) ,DELN(600) ,DEL(600) ,GNM1 (600) ,GNM2(600) 
C0MM0N/MASS/A(600) ,B(600) ,CM(8) 

COMMON /DELT/XDEL ( 4 , 1 8 ) 

COMMON /F0R/FF(5) ,FC(5) 

COMMON / SIG  Z /DSIG  Z ( 1 2 ) 

COMMON /FIRST /F0( 600 ) ,DER( 600) 

DIMENSION  AEI(12) 

IF(NTIME.NE. 1 )G0  TO  15 
DO  16  1  =  1  ,6 
16  DSIGZ(I  )=0 .0 
GO  TO  110 
15  CONTINUE 
AEIT=0 . 

DO  20  1=1,4 
XXX(I)=X(IX(N ,1) ) 

YYY(I)=Y(IX(N,I)) 

20  CONTINUE 

XXX(5)=(XXX(1 )+XXX(2)+XXX(3)+XXX(4))/4.0 
YYY(5)=(YYY(1 )+YYY(2)+YYY(3 )+YYY(4) )/4 .0 
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FAV=0 . 

CALL  DISFOR 
DO  30  1=1,5 
30  FAV=FAV+FF(I)/5. 

DO  100  11=1 , NUMLA 

C  COMPUTE  AVERAGE  ELEMENT  ACCELERATION  (AEA) 

I=NUMLA+1-II 

AEC1=0. 

DO  40  J=1 ,4 
DO  40  K=1 ,2 

INL=IX(N,J)*3+NUMNP*3*(I+1-K) 

IF(NTIME.EQ. 1 )GO  TO  35 

AEC1=AEC1+(DEL(INL)-2.*DELN(INL)+DELN1(INL))/(8.*H**2) 
GO  TO  40 

35  AEC1=AEC1+((DEL(INL)-DELN(INL))/H-DER(INL))/(8.*H) 

40  CONTINUE 
AEI(I )=AEC1 

100  AEIT=AEIT+AEC1 /FLOAT (NUMLA) 

C  COMPUTE  DIFFERENCE  BETWEEN  INDIVIDUAL  ELEMENT  ACCELERATION 
C  AND  TOTAL  PLATE  ELEMENT  ACCELERATION 
SART=0 . 

DO  50  1=1, NUMLA 
MTR=MATRIL(I) 

RT=TH(I)*RO(MTR) 

SART=SART+RT*AEI (I) 

50  CONTINUE 

IF ( SART . EQ . 0 . 0 )GO  TO  52 

AK=FAV/SART 

GO  TO  53 

52  AK=1 . 

53  CONTINUE 

C  SOLVE  FOR  SIGMAZ 
SN=FAV 

DO  80  11=1, NUMLA 
I=NUMLA+1-II 
MTR=MATRIL(I ) 

SNM=SN-AK«AEI(I)«RO(MTR)*TH(I) 

SIGZ= (SNM+SN)/2 .0 

DSIGZ ( I ) =SIGZ-SIGMA (I ,N , 3) 

80  SN=SNM 
110  CONTINUE 
RETURN 
END 


j 
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VISCOPLASTIC  YIELD  SUBROUTINE 
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SUBROUTINE  YIELD 
COMPLEX  LAMBDA,  DISC,  DISQ 

COMMON /ARG/XXX( 10) ,YYY(10) ,S(24) ,XX(3) ,YY(3) , 

1 CRZ  ( 6 , 6 )  ,  XI  ( 1 0 ) ,  SIG  ( 1 2 ) ,  N ,  M 
C0MM0N/3IGM/BSUM,TSUM(6) , SIGMA (12,50,6) , DSIG (6) ,F(600) 

COMMON /BASIC2/BET , H , HN , HH ( 3 ) , HT ( 2 ) , TIME , NLAY , IF LAG , IT , NTIME 
C0MM0N/PLYLD/SIGY(12,6) ,DEPS(6) 

COMMON / MAT P/RO( 12) ,E(9 , 12) ,EE(9) ,ETA( 12) 

COMMON / SIGZ/DSIGZ ( 1 2 ) 

DIMENSION  SIGT(6)  ,SIGYB(3)  ,TSIG(6) ,TSIGB(6) 
SIG3=SIGMA(M,N,3)+DSIGZ(M) 

C  FIND  COMBINATION  YIELD  STRESSES  (SIGYB) 

SIGYB( 1 ) = 1 . O/SIGY (M, 1 )**2-1 .0/SIGY(M,2)**2-1 . O/SIGY (M, 3 )**2 
SIGYB(2)=1 . O/SIGY (M, 2 )**2-1 . O/SIGY (M, 1 ) **2-1 . O/SIGY (M, 3 )**2 
SIGYB(3)=1 .O/SIGY (M,3)**2-1 .0/SIGY(M, 1 )**2-1 . O/SIGY (M,2)**2 
C  GET  TEST  STRESSES  FROM  ELASTIC  ANALYSIS  (SIGT) 

DO  15  1=1,6 

15  SIGT (I) =SIGMA(M,N , I)+DSIG(I ) 

C  TEST  FOR  YIELDING  BY  HILL'S  CRITERION 
YLD=0 .0 
DO  20  1  =  1  ,6 

20  YLD=YLD+SIGT(I)**2/3IGY(M,I)**2 

YLD=YLD+SIG YB ( 1 ) *  SIGT ( 2 ) *SIGT ( 3 ) +SIGYB ( 2 ) *SIGT ( 1 ) *3IGT ( 3 ) 

1  +SIGYB(3)*SIGT(1 )*SIGT(2) 

IF(YLD.LT. 1 .O.OR.YLD.EQ. 1.0)  GO  TO  390 
C  CORRECT  STRESSES  FUR  PLASTICITY 

C  FIND  TEST  STRAIN  INCREMENT  VALUES  (TSIG  *  LAMBDA  =  DEP3  PLASTIC) 
TSIG ( 1 ) =SIGT ( 1 ) /SIGY (M, 1 )**2+(SIGYB( 3 )*SIGT (2 )+SlGYB( 2 )* 

1  SIGT(3))/2.0 

TSIG (2 ) =SIGT ( 2 ) /SIGY (M , 2 ) *  *2+ (SIGYB ( 3 ) *SIGT ( 1 ) +SIGYB ( 1 ) * 

1  SIGT (3 ) ) /2 .0 

TSIG(3)=SIGT(3)/SIGY(M,3)**2+(S1GYB(2)*SIGT(1 )+SIGYB(1 )* 

1  SIGT (2 ) ) /2 .0 
DO  25  1=4,6 

25  TSIG(I)=SIGT(I)/(SIGY (M,I)**2) 

DO  35  1=1,6 
TSIGB(I)=0.0 
DO  35  J=1 ,6 

35  TSIGB(I )=TSIGB(I )+CRZ(I , J )*TSIG( J ) 

C  FIND  A ,B, C  OF  THE  QUADRADIC  EQUATION  TO  CALCULATE  LAMBDA 
AQ=0 .0 
DO  45  1=1,6 

45  AQ=AQ+(TSIGB(I )+ETA(M)*SIGT(I)/H)**2/SIGY(M,I)**2 

AQ=AQ+SIGYB(1 )*(TSIGB(2)+ETA(M)*SIGT(2)/H)*(TSIGB(3)+ETA(M) 

1  *SIGT(3)/H)+SIGYB(2)*(TSIGB(1 )+ETA(M)*SIGT(1 )/H)*(TSIGB(3) 

2  +ETA(M)*SIGT(3)/H)+SIGYB(3)*(TSIGB(1 )+ETA(M)*SIGT(1)/H)* 

3  (TSIGB ( 2 ) +ETA ( M) *3IGT ( 2 ) /H ) 

BQ=0.0 

DO  55  1  =  1  ,6 

55  BQ=BQ+(TSIGB(I)+ETA(M)*SIGT(I )/H)*SIGT(I )/SIGY(M, I)**2 
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BQ=BQ+ (SIGYB ( 1 ) • ( SIGT (2 ) • ( TSIGB ( 3 ) +ETA (M)*SIGT ( 3 ) /H ) +SIGT ( 3 ) 

1  «(TSIGB(2)+ETA(M)»SIGT(2)/H))+SIGYB(2)»(SIGT(1)«(TSIGB(3) 

2  +ETA ( M) •SIGT (3)/H) +SIGT (3) *( TSIGB ( 1 ) +ETA (M) *SIGT (1 )/H))+ 

3  SIGYB ( 3 ) • ( SIGT ( 1 ) • ( TSIGB ( 2 ) +ETA ( M) »SIGT ( 2 ) /H ) +SIGT ( 2 ) • ( TSIGB ( 1 ) 

4  +ETA(M)»SIGTCl)/H)))/2.0 
CQ=YLD-1 .0 
DISC=BQ««2-AQ»CQ 
DISQ=CSQRT (DISC) 

LAMBDA=CQ/ ( BQ+DISQ ) 

C  CALCULATE  ELASTIC-PLASTIC  STRESSES 
DO  65  1=1,6 

65  SIGMA(M,N,I)=(1 . 0+L AMBDA*ETA ( M) /H ) • ( SIGT ( I ) -LAMBDA* ( TSIGB ( I ) 

1  +ETA(M)*SIGT(I)/H) ) 

GO  TO  400 
390  DO  395  1=1,6 
395  SIGMA ( M , N , I ) =SIGT ( I ) 

GO  TO  600 
400  CONTINUE 

SIG3 = ( SIG  3+SIGMA ( M , N , 3 ) ) /2 . 

SIGMA (M,N,2)=SIGMA(M,N,2)-SIGMA(M,N,3)+SIG3 
SIGMA ( M , N , 1 ) =SIGMA (M,N , 1 ) -SIGMA ( M , N , 3 ) +SIG  3 
SIGMA ( M , N , 3 ) =SIG  3 
600  CONTINUE 
RETURN 
END 
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EXAMPLE  DATA  INPUT 
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EXAMPLE - STEEL-ALUMINUM-STEEL  18X9  INCHES  05/29/79  DAT  1C 

3  20  1 


2 

1 


.25 

6 


1 

1 

6 

6 


2.50E-06 


0.0 

4.5 

4.5 

0.0 


1  1 

2  3 

4  4 

5  5 

8  8 

9  9 

12  12 

13  13 

16  16 

17  17 

20  20 

21  21 

22  23 

24  24 

25  25 

26  28 

29  29 

33  33 

37  37 

41  41 

45  45 

49  49 

50  52 

53  53 

57  57 

61  61 

65  65 

69  69 

73  73 

74  76 

77  77 

81  81 

85  85 

89  89 

93  93 

1 

•3500000E+05 

2 

.5000000E+05 

1 

.3500000E+05 


35  v 
0.0 
0.0 
9.0 
9.0 

4 
2 

5 
1 
3 
1 
3 
1 
3 
1 
3 

6 

3 
7 

4 
2 
1 
1 
1 
1 
1 
4 
2 
1 
1 
1 
1 
1 
4 
2 
1 
1 
1 
1 
1 


5.00E-06 

3 

4  1 

4  6 

1  6 
1  1 


1.00E-05  2.00E-05  6.00E-05 


4.5 

4.5 

0.0 

0.0 


0.0 

9.0 

9.0 

0.0 


6.243388E+00 
3.913917  E— 0 1 
6.243388E+00 
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1  .00073236 

3.0E+07  3.0E+07  3.0E+07  .25  .25  .25 

1.2E+07  1.2E+07  1.2E+07 

2  .00025251 

1.0E+07  1.0E+07’  1.0E+07  .3  .3  .3 

.3846E+07  .3846E+07  .3846E+07 
.25  .5  .5 

GRID  SIZE  FOR  THIS  CASE  IS  4  X  6  ,  4/29/79,  DAT 1C 
4.5  9.0 
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